model/src/ini_depths.F

C $Header: /u/gcmpack/MITgcm/model/src/ini_depths.F,v 1.26 2001/09/26 18:09:15 cnh 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  - Loop counters
C     I,J,K
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 (groundAtK1 .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 (groundAtK1 .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
        _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)

      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 ( groundAtK1 .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( topoZ, Ro_surf, myThid )

          _BARRIER
          _BEGIN_MASTER( myThid )
          CALL WRITE_FLD_XY_RS( 'topo_P',' ',Ro_surf,0,myThid) 
          _END_MASTER(myThid)

        ELSE
C----
C   Direct Transfer to Ro_surf :
          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 (groundAtK1) 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	jmc	1.27	C $Header: /u/gcmpack/MITgcm/model/src/ini_depths.F,v 1.26 2001/09/26 18:09:15 cnh 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			C \| SUBROUTINE INI_DEPTHS
13			C \| o define R_position of Lower and Surface Boundaries
14			C ==========================================================
15			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			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			C bi,bj - Loop counters
50	adcroft	1.21	C I,J,K
51	cnh	1.13	C oldPrec - Temporary used in controlling binary input dataset precision
52	jmc	1.24	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			IF (groundAtK1 .AND. bathyFile .NE. ' '
60			& .AND. topoFile .NE. ' ' ) THEN
61			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.25	IF (groundAtK1 .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	jmc	1.24	_BEGIN_MASTER( myThid )
100	adcroft	1.19	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.
102	jmc	1.24	CALL READ_REC_XY_RS( bathyFile, R_low, 1, 0, myThid )
103	adcroft	1.19	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
105	jmc	1.24	c CALL READ_FLD_XY_RS( bathyFile, ' ', R_low, 0, myThid )
106	adcroft	1.19	C Read the bathymetry using the low-level I/O package
107	jmc	1.24	c CALL MDSREADFIELD( bathyFile, readBinaryPrec,
108			c & 'RS', 1, R_low, 1, myThid )
109			_END_MASTER(myThid)
110
111	jmc	1.25	ENDIF
112			C- end setup R_low in the interior
113	jmc	1.24
114	jmc	1.25	C- fill in the overlap :
115			_EXCH_XY_R4(R_low, myThid )
116	adcroft	1.21
117	jmc	1.24	c CALL PLOT_FIELD_XYRS(R_low,'Bottom depths (ini_depths)',1,myThid)
118			c _BEGIN_MASTER( myThid )
119			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	jmc	1.25	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	jmc	1.24
139			DO bj = myByLo(myThid), myByHi(myThid)
140			DO bi = myBxLo(myThid), myBxHi(myThid)
141	jmc	1.25	DO j=1,sNy
142			DO i=1,sNx
143	jmc	1.24	Ro_surf(i,j,bi,bj) = Ro_SeaLevel
144			ENDDO
145			ENDDO
146			ENDDO
147			ENDDO
148
149			ELSE
150	jmc	1.25	C------ read from file:
151	jmc	1.24
152			C- read surface topography (in m) from topoFile (case topoFile.NE.' '):
153			_BEGIN_MASTER( myThid )
154	jmc	1.27	CALL READ_REC_XY_RS( topoFile, topoZ, 1, 0, myThid )
155	jmc	1.24	_END_MASTER(myThid)
156			_BARRIER
157
158			IF (buoyancyRelation .EQ. 'ATMOSPHERIC') THEN
159	jmc	1.25	C----
160	jmc	1.24	C Convert Surface Geopotential to (reference) Surface Pressure
161			C according to Tref profile, using same discretisation as in calc_phi_hyd
162	jmc	1.25	C----
163	jmc	1.24	c _BEGIN_MASTER( myThid )
164	jmc	1.27	c CALL WRITE_FLD_XY_RS( 'topo_Z',' ',topoZ,0,myThid)
165	jmc	1.24	c _END_MASTER(myThid)
166
167	jmc	1.27	CALL INI_P_GROUND( topoZ, Ro_surf, myThid )
168	jmc	1.24
169			_BARRIER
170			_BEGIN_MASTER( myThid )
171			CALL WRITE_FLD_XY_RS( 'topo_P',' ',Ro_surf,0,myThid)
172			_END_MASTER(myThid)
173
174			ELSE
175	jmc	1.25	C----
176	jmc	1.24	C Direct Transfer to Ro_surf :
177			DO bj = myByLo(myThid), myByHi(myThid)
178			DO bi = myBxLo(myThid), myBxHi(myThid)
179			DO j=1,sNy
180			DO i=1,sNx
181	jmc	1.27	Ro_surf(i,j,bi,bj) = topoZ(i,j,bi,bj)
182	jmc	1.24	ENDDO
183			ENDDO
184			ENDDO
185			ENDDO
186	jmc	1.25
187	jmc	1.24	ENDIF
188
189	jmc	1.25	C------ end case "read topoFile"
190			ENDIF
191	jmc	1.24
192	jmc	1.25	C----- fill in the overlap :
193			_EXCH_XY_R4(Ro_surf, myThid )
194	jmc	1.24
195			c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
196	adcroft	1.21
197	jmc	1.24	C------
198			C 3) Close the Domain (special configuration).
199			C------
200	adcroft	1.23	IF (groundAtK1) THEN
201	adcroft	1.21	DO bj = myByLo(myThid), myByHi(myThid)
202			DO bi = myBxLo(myThid), myBxHi(myThid)
203			DO j=1-Oly,sNy+Oly
204			DO i=1-Olx,sNx+Olx
205	jmc	1.24	iG = myXGlobalLo-1+(bi-1)*sNx+I
206			jG = myYGlobalLo-1+(bj-1)*sNy+J
207			C Test for eastern edge
208			c IF ( iG .EQ. Nx ) Ro_surf(i,j,bi,bj) = 0.
209			C Test for northern edge
210			c IF ( jG .EQ. Ny ) Ro_surf(i,j,bi,bj) = 0.
211	adcroft	1.23	IF (usingSphericalPolarGrid .AND. abs(yC(I,J,bi,bj)).GE.90. )
212	jmc	1.24	& Ro_surf(I,J,bi,bj) = rF(Nr+1)
213	adcroft	1.23	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.24	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.23	IF (usingSphericalPolarGrid .AND. abs(yC(I,J,bi,bj)).GE.90. )
229	jmc	1.24	& R_low(I,J,bi,bj) = Ro_SeaLevel
230	adcroft	1.21	ENDDO
231			ENDDO
232			ENDDO
233			ENDDO
234			ENDIF
235	jmc	1.24
236			c _BEGIN_MASTER( myThid )
237			c CALL WRITE_FLD_XY_RS('Ro_surf',' ',Ro_surf,0,myThid)
238			c _END_MASTER(myThid)
239	adcroft	1.23
240	cnh	1.1	RETURN
241			END