model/src/ini_depths.F

C $Header: /u/u0/gcmpack/MITgcm/model/src/ini_depths.F,v 1.31 2005/06/22 00:18:00 jmc Exp $
C $Name:  $

#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: INI_DEPTHS
C     !INTERFACE:
      SUBROUTINE INI_DEPTHS( myThid )
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE INI_DEPTHS
C     | o define R_position of Lower and Surface Boundaries
C     *==========================================================*
C     |atmosphere orography:
C     | define either in term of P_topo or converted from Z_topo
C     |ocean bathymetry:
C     | The depths of the bottom of the model is specified in
C     | terms of an XY map with one depth for each column of
C     | grid cells. Depths do not have to coincide with the
C     | model levels. The model lopping algorithm makes it
C     | possible to represent arbitrary depths.
C     | The mode depths map also influences the models topology
C     | By default the model domain wraps around in X and Y.
C     | This default doubly periodic topology is "supressed"
C     | if a depth map is defined which closes off all wrap
C     | around flow.
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SURFACE.h"

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     myThid -  Number of this instance of INI_DEPTHS
      INTEGER myThid
CEndOfInterface

C     !LOCAL VARIABLES:
C     == Local variables ==
C     iG, jG - Global coordinate index
C     bi,bj  - Tile indices
C     I,J,K  - Loop counters
C     oldPrec - Temporary used in controlling binary input dataset precision
C     msgBuf    - Informational/error meesage buffer
      INTEGER iG, jG
      INTEGER bi, bj
      INTEGER  I, J
      CHARACTER*(MAX_LEN_MBUF) msgBuf
CEOP

      IF (usingPCoords .AND. bathyFile .NE. ' '
     &                 .AND. topoFile  .NE. ' ' ) THEN
       WRITE(msgBuf,'(A,A)')
     &  'S/R INI_DEPTHS: both bathyFile & topoFile are specified:',
     &  ' select the right one !'
       CALL PRINT_ERROR( msgBuf , myThid)
       STOP 'ABNORMAL END: S/R INI_DEPTHS'
      ENDIF

C------
C   0) Initialize R_low and Ro_surf (define an empty domain)
C------
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO j=1-Oly,sNy+Oly
         DO i=1-Olx,sNx+Olx
          R_low(i,j,bi,bj) = 0.
          Ro_surf(i,j,bi,bj) = 0.
          topoZ(i,j,bi,bj) = 0.
         ENDDO
        ENDDO
       ENDDO
      ENDDO

C------
C   1) Set R_low = the Lower (in r sense) boundary of the fluid column :
C------
      IF (usingPCoords .OR. bathyFile .EQ. ' ') THEN
C- e.g., atmosphere : R_low = Top of atmosphere
C-            ocean : R_low = Bottom
       DO bj = myByLo(myThid), myByHi(myThid)
        DO bi = myBxLo(myThid), myBxHi(myThid)
         DO j=1,sNy
          DO i=1,sNx
           R_low(i,j,bi,bj) = rF(Nr+1)
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ELSE
        _BARRIER
        _BEGIN_MASTER( myThid )
C Read the bathymetry using the mid-level I/O pacakage read_write_rec
C The 0 is the "iteration" argument. The 1 is the record number.
        CALL READ_REC_XY_RS( bathyFile, R_low, 1, 0, myThid )
C Read the bathymetry using the mid-level I/O pacakage read_write_fld
C The 0 is the "iteration" argument. The ' ' is an empty suffix
c       CALL READ_FLD_XY_RS( bathyFile, ' ', R_low, 0, myThid )
C Read the bathymetry using the low-level I/O package
c       CALL MDSREADFIELD( bathyFile, readBinaryPrec,
c    &                     'RS', 1, R_low, 1, myThid )
        _END_MASTER(myThid)
        _BARRIER
      ENDIF
C- end setup R_low in the interior

C- fill in the overlap :
      _EXCH_XY_R4(R_low, myThid )

c     CALL PLOT_FIELD_XYRS(R_low,'Bottom depths (ini_depths)',1,myThid)
c     _BEGIN_MASTER( myThid )
c     CALL WRITE_FLD_XY_RS( 'R_low' ,' ', R_low, 0,myThid)
c     _END_MASTER(myThid)

c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C------
C   2) Set R_surf = Surface boundary: ocean surface / ground for the atmosphere
C------

      IF ( usingPCoords .AND. bathyFile.NE.' ' ) THEN
C------ read directly Po_surf from bathyFile (only for backward compatibility)

        _BEGIN_MASTER( myThid )
        CALL READ_REC_XY_RS( bathyFile, Ro_surf, 1, 0, myThid )
        _END_MASTER(myThid)
        _BARRIER

      ELSEIF ( topoFile.EQ.' ' ) THEN
C------ set default value:

        DO bj = myByLo(myThid), myByHi(myThid)
         DO bi = myBxLo(myThid), myBxHi(myThid)
          DO j=1,sNy
           DO i=1,sNx
            Ro_surf(i,j,bi,bj) = Ro_SeaLevel
           ENDDO
          ENDDO
         ENDDO
        ENDDO

      ELSE
C------ read from file:

C- read surface topography (in m) from topoFile (case topoFile.NE.' '):
        _BEGIN_MASTER( myThid )
        CALL READ_REC_XY_RS( topoFile, topoZ, 1, 0, myThid )
        _END_MASTER(myThid)
        _BARRIER

        IF (buoyancyRelation .EQ. 'ATMOSPHERIC') THEN
C----
C   Convert Surface Geopotential to (reference) Surface Pressure
C   according to Tref profile, using same discretisation as in calc_phi_hyd
C----
c         _BEGIN_MASTER( myThid )
c         CALL WRITE_FLD_XY_RS( 'topo_Z',' ',topoZ,0,myThid)
c         _END_MASTER(myThid)

          CALL INI_P_GROUND( 2, topoZ,
     O                       Ro_surf,
     I                       myThid )

          _BARRIER
C         This I/O is now done in write_grid.F
c         _BEGIN_MASTER( myThid )
c         CALL WRITE_FLD_XY_RS( 'topo_P',' ',Ro_surf,0,myThid)
c         _END_MASTER(myThid)

        ELSE
C----
C   Direct Transfer to Ro_surf (e.g., to specify upper ocean boundary
C    below an ice-shelf - NOTE - actually not yet implemented )
          DO bj = myByLo(myThid), myByHi(myThid)
           DO bi = myBxLo(myThid), myBxHi(myThid)
            DO j=1,sNy
             DO i=1,sNx
              Ro_surf(i,j,bi,bj) = topoZ(i,j,bi,bj)
             ENDDO
            ENDDO
           ENDDO
          ENDDO

        ENDIF

C------ end case "read topoFile"
      ENDIF

C----- fill in the overlap :
      _EXCH_XY_R4(Ro_surf, myThid )

c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C------
C   3) Close the Domain (special configuration).
C------
      IF (usingPCoords) THEN
       DO bj = myByLo(myThid), myByHi(myThid)
        DO bi = myBxLo(myThid), myBxHi(myThid)
         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx
           iG = myXGlobalLo-1+(bi-1)*sNx+I
           jG = myYGlobalLo-1+(bj-1)*sNy+J
C          Test for eastern edge
c          IF ( iG .EQ. Nx )  Ro_surf(i,j,bi,bj) = 0.
C          Test for northern edge
c          IF ( jG .EQ. Ny )  Ro_surf(i,j,bi,bj) = 0.
           IF (usingSphericalPolarGrid .AND. abs(yC(I,J,bi,bj)).GE.90. )
     &       Ro_surf(I,J,bi,bj) = rF(Nr+1)
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ELSE
       DO bj = myByLo(myThid), myByHi(myThid)
        DO bi = myBxLo(myThid), myBxHi(myThid)
         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx
           iG = myXGlobalLo-1+(bi-1)*sNx+I
           jG = myYGlobalLo-1+(bj-1)*sNy+J
C          Test for eastern edge
c          IF ( iG .EQ. Nx )  R_low(i,j,bi,bj) = 0.
C          Test for northern edge
c          IF ( jG .EQ. Ny )  R_low(i,j,bi,bj) = 0.
           IF (usingSphericalPolarGrid .AND. abs(yC(I,J,bi,bj)).GE.90. )
     &       R_low(I,J,bi,bj) = Ro_SeaLevel
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDIF

c     _BEGIN_MASTER( myThid )
c     CALL WRITE_FLD_XY_RS('Ro_surf',' ',Ro_surf,0,myThid)
c     _END_MASTER(myThid)

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