aim.5l_Equatorial_Channel/code/ini_depths.F

C $Header: /u/gcmpack/models/MITgcmUV/verification/aim.5l_Equatorial_Channel/code/ini_depths.F,v 1.3.2.2 2001/05/04 13:27:04 jmc Exp $
C $Name:  $

#include "CPP_OPTIONS.h"

CStartOfInterface
      SUBROUTINE INI_DEPTHS( myThid )
C     /==========================================================\
C     | SUBROUTINE INI_DEPTHS                                    |
C     | o Initialise map of model depths                         |
C     |==========================================================|
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     \==========================================================/
      IMPLICIT NONE

C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "INI_DEPTHS.h"

C     == Routine arguments ==
C     myThid -  Number of this instance of INI_DEPTHS
      INTEGER myThid
CEndOfInterface

C     == Local variables ==
C     iG, jG - Global coordinate index
C     bi,bj  - Loop counters
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--------------------
C     oldPrec - Temporary used in controlling binary input dataset precision
      INTEGER iG, jG
      INTEGER bi, bj
      INTEGER  I, J
      _RL Hdefault

      _BARRIER
      IF ( bathyFile .EQ. ' '  ) THEN
C      Set up a flat bottom box with doubly periodic topology.
C      H is the basic variable from which other terms are derived. It
C      is the term that would be set from an external file for a 
C      realistic problem.
       IF (groundAtK1) THEN
         Hdefault = Ro_SeaLevel
       ELSE
         Hdefault = rF(Nr+1)
       ENDIF
       DO bj = myByLo(myThid), myByHi(myThid)
        DO bi = myBxLo(myThid), myBxHi(myThid)
         DO j=1,sNy
          DO i=1,sNx
           iG = myXGlobalLo-1+(bi-1)*sNx+I
           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.
          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, H, 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, ' ', H, 0, myThid )
C Read the bathymetry using the low-level I/O package
C      CALL MDSREADFIELD( bathyFile, readBinaryPrec,
C    &                    'RS', 1, H, 1, myThid )
       _END_MASTER(myThid)
      ENDIF

      _EXCH_XY_R4(    H, myThid )

C
C     CALL PLOT_FIELD_XYRS(H,'Model depths (ini_depths)',1,myThid)
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
           R_low(I,J,bi,bj) = rF(Nr+1)
           Ro_surf(I,J,bi,bj) = H(I,J,bi,bj)
c- Domain : Symetric % Eq. & closed at N & S boundaries:
           IF (usingSphericalPolarGrid .AND.
     &         abs(yC(I,J,bi,bj)).GE.-phiMin) Ro_surf(I,J,bi,bj) = 0.
           IF (usingSphericalPolarGrid .AND. abs(yC(I,J,bi,bj)).GE.90. ) 
     &     THEN
            Ro_surf(I,J,bi,bj) = 0.
           ENDIF
          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
           R_low(I,J,bi,bj) = H(I,J,bi,bj)
           Ro_surf(I,J,bi,bj) = Ro_SeaLevel
c- Domain : Symetric % Eq. & closed at N & S boundaries:
           IF (usingSphericalPolarGrid .AND.
     &         abs(yC(I,J,bi,bj)).GE.-phiMin) R_low(I,J,bi,bj) = 0.
           IF (usingSphericalPolarGrid .AND. abs(yC(I,J,bi,bj)).GE.90. ) 
     &     THEN
            R_low(I,J,bi,bj) = 0.
           ENDIF
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDIF

      RETURN
      END
1	C $Header: /u/gcmpack/models/MITgcmUV/verification/aim.5l_Equatorial_Channel/code/ini_depths.F,v 1.3.2.2 2001/05/04 13:27:04 jmc Exp $
2	C $Name: $
3
4	#include "CPP_OPTIONS.h"
5
6	CStartOfInterface
7	SUBROUTINE INI_DEPTHS( myThid )
8	C /==========================================================\
9	C \| SUBROUTINE INI_DEPTHS \|
10	C \| o Initialise map of model depths \|
11	C \|==========================================================\|
12	C \| The depths of the bottom of the model is specified in \|
13	C \| terms of an XY map with one depth for each column of \|
14	C \| grid cells. Depths do not have to coincide with the \|
15	C \| model levels. The model lopping algorithm makes it \|
16	C \| possible to represent arbitrary depths. \|
17	C \| The mode depths map also influences the models topology \|
18	C \| By default the model domain wraps around in X and Y. \|
19	C \| This default doubly periodic topology is "supressed" \|
20	C \| if a depth map is defined which closes off all wrap \|
21	C \| around flow. \|
22	C \==========================================================/
23	IMPLICIT NONE
24
25	C === Global variables ===
26	#include "SIZE.h"
27	#include "EEPARAMS.h"
28	#include "PARAMS.h"
29	#include "GRID.h"
30	#include "INI_DEPTHS.h"
31
32	C == Routine arguments ==
33	C myThid - Number of this instance of INI_DEPTHS
34	INTEGER myThid
35	CEndOfInterface
36
37	C == Local variables ==
38	C iG, jG - Global coordinate index
39	C bi,bj - Loop counters
40	C I,J,K
41	C H - Depth of base of fluid from upper surface f[X,Y] (m).
42	C Hdefault - default r-coordinate of the lower boundary (=ground)
43	C (=minus(Total_depth) in the ocean model)
44	C (=Total Pressure at Sea Level in the atmos model)
45	C--------------------
46	C NOTE: will change soon: 2 separed files for r_lower and r_surface boudaries
47	C and for the atmosphere, topography will be defined in term of height
48	C--------------------
49	C oldPrec - Temporary used in controlling binary input dataset precision
50	INTEGER iG, jG
51	INTEGER bi, bj
52	INTEGER I, J
53	_RL Hdefault
54
55	_BARRIER
56	IF ( bathyFile .EQ. ' ' ) THEN
57	C Set up a flat bottom box with doubly periodic topology.
58	C H is the basic variable from which other terms are derived. It
59	C is the term that would be set from an external file for a
60	C realistic problem.
61	IF (groundAtK1) THEN
62	Hdefault = Ro_SeaLevel
63	ELSE
64	Hdefault = rF(Nr+1)
65	ENDIF
66	DO bj = myByLo(myThid), myByHi(myThid)
67	DO bi = myBxLo(myThid), myBxHi(myThid)
68	DO j=1,sNy
69	DO i=1,sNx
70	iG = myXGlobalLo-1+(bi-1)*sNx+I
71	jG = myYGlobalLo-1+(bj-1)*sNy+J
72	C Default depth of full domain
73	H(i,j,bi,bj) = Hdefault
74	C Test for eastern edge
75	C IF ( iG .EQ. nX ) H(i,j,bi,bj) = 0.
76	C Test for northern edge
77	C IF ( jG .EQ. nY ) H(i,j,bi,bj) = 0.
78	ENDDO
79	ENDDO
80	ENDDO
81	ENDDO
82	ELSE
83	_BEGIN_MASTER( myThid )
84	C Read the bathymetry using the mid-level I/O pacakage read_write_rec
85	C The 0 is the "iteration" argument. The 1 is the record number.
86	CALL READ_REC_XY_RS( bathyFile, H, 1, 0, myThid )
87	C Read the bathymetry using the mid-level I/O pacakage read_write_fld
88	C The 0 is the "iteration" argument. The ' ' is an empty suffix
89	C CALL READ_FLD_XY_RS( bathyFile, ' ', H, 0, myThid )
90	C Read the bathymetry using the low-level I/O package
91	C CALL MDSREADFIELD( bathyFile, readBinaryPrec,
92	C & 'RS', 1, H, 1, myThid )
93	_END_MASTER(myThid)
94	ENDIF
95
96	_EXCH_XY_R4( H, myThid )
97
98	C
99	C CALL PLOT_FIELD_XYRS(H,'Model depths (ini_depths)',1,myThid)
100	C
101	IF (groundAtK1) THEN
102	DO bj = myByLo(myThid), myByHi(myThid)
103	DO bi = myBxLo(myThid), myBxHi(myThid)
104	DO j=1-Oly,sNy+Oly
105	DO i=1-Olx,sNx+Olx
106	R_low(I,J,bi,bj) = rF(Nr+1)
107	Ro_surf(I,J,bi,bj) = H(I,J,bi,bj)
108	c- Domain : Symetric % Eq. & closed at N & S boundaries:
109	IF (usingSphericalPolarGrid .AND.
110	& abs(yC(I,J,bi,bj)).GE.-phiMin) Ro_surf(I,J,bi,bj) = 0.
111	IF (usingSphericalPolarGrid .AND. abs(yC(I,J,bi,bj)).GE.90. )
112	& THEN
113	Ro_surf(I,J,bi,bj) = 0.
114	ENDIF
115	ENDDO
116	ENDDO
117	ENDDO
118	ENDDO
119	ELSE
120	DO bj = myByLo(myThid), myByHi(myThid)
121	DO bi = myBxLo(myThid), myBxHi(myThid)
122	DO j=1-Oly,sNy+Oly
123	DO i=1-Olx,sNx+Olx
124	R_low(I,J,bi,bj) = H(I,J,bi,bj)
125	Ro_surf(I,J,bi,bj) = Ro_SeaLevel
126	c- Domain : Symetric % Eq. & closed at N & S boundaries:
127	IF (usingSphericalPolarGrid .AND.
128	& abs(yC(I,J,bi,bj)).GE.-phiMin) R_low(I,J,bi,bj) = 0.
129	IF (usingSphericalPolarGrid .AND. abs(yC(I,J,bi,bj)).GE.90. )
130	& THEN
131	R_low(I,J,bi,bj) = 0.
132	ENDIF
133	ENDDO
134	ENDDO
135	ENDDO
136	ENDDO
137	ENDIF
138
139	RETURN
140	END