model/src/ini_depths.F

C $Header: /u/gcmpack/models/MITgcmUV/model/src/ini_depths.F,v 1.2 1998/04/24 02:05:41 cnh Exp $

#include "CPP_EEOPTIONS.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's 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     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"

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

C     == Local variables ==
C     xG, yG - Global coordinate location.
C     zG
C     zUpper - Work arrays for upper and lower 
C     zLower   cell-face heights.
C     phi    - Temporary scalar
C     iG, jG - Global coordinate index
C     bi,bj  - Loop counters
C     zUpper - Temporary arrays holding z coordinates of
C     zLower   upper and lower faces.
C     I,J,K
      _RL    xG, yG, zG
      _RL    phi
      _RL    zUpper(Nz), zLower(Nz)
      INTEGER iG, jG
      INTEGER bi, bj
      INTEGER  I,  J, K

C     For now 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.
      _BARRIER
      phi = 0. D0
      DO K=1,Nz
       phi = phi+delZ(K)
      ENDDO
      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) = phi
C         Test for eastern edge
          IF ( iG .EQ. nX ) H(i,j,bi,bj) = 0.
C         Test for northern edge
          IF ( jG .EQ. nY ) H(i,j,bi,bj) = 0.
C         Island
          IF ( iG .EQ. 1 .AND. 
     &         jG .EQ. 24 ) H(i,j,bi,bj) = 0.75*phi
C         IF ( iG .EQ. 1 .AND. 
C    &         jG .EQ. 24 ) H(i,j,bi,bj) = 0.
           
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO J=1,sNy
         DO I=1,sNx
C         Inverse of depth
          IF ( h(i,j,bi,bj) .EQ. 0. _d 0 ) THEN
           rH(i,j,bi,bj) = 0. _d 0
          ELSE
           rH(i,j,bi,bj) = 1. _d 0 /  h(i,j,bi,bj)
          ENDIF
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      _EXCH_XY_R4(    H, myThid )
      _EXCH_XY_R4(   rH, myThid )

C--   Now calculate "lopping" factors hFac.
      zG = delZ(1)*0.5 D0
      zFace(1) = 0
      rDzC(1)  = 2. _d 0/delZ(1)
      DO K=1,Nz-1
       saFac(K) = 1. D0
       zC(K)      = zG
       zG         = zG + (delZ(K)+delZ(K+1))*0.5 D0
       zFace(K+1) = zFace(K)+delZ(K)
       rDzC(K+1)  = 2. _d 0/(delZ(K)+delZ(K+1))
      ENDDO
      zC(Nz)      = zG
      zFace(Nz+1) = zFace(Nz)+delZ(Nz)
      DO K=1,Nz
       zUpper(K) = zFace(K)
       zLower(K) = zFace(K+1)
       dzF(K)    = delz(K)
       rdzF(K)   = 1.d0/dzF(K)
      ENDDO
      DO bj=myByLo(myThid), myByHi(myThid)
       DO bi=myBxLo(myThid), myBxHi(myThid)
        DO K=1, Nz
         DO J=1,sNy
          DO I=1,sNx
           hFacC(I,J,K,bi,bj) = 1.
           IF     ( H(I,J,bi,bj) .LE. zUpper(K) ) THEN
C           Below base of domain
            hFacC(I,J,K,bi,bj) = 0.
           ELSEIF ( H(I,J,bi,bj) .GT. zLower(K) ) THEN
C           Base of domain is below this cell
            hFacC(I,J,K,bi,bj) = 1.
           ELSE
C           Base of domain is in this cell
C           Set hFac tp the fraction of the cell that is open.
            phi = zUpper(K)-H(I,J,bi,bj)
            hFacC(I,J,K,bi,bj) = phi/(zUpper(K)-zLower(K))
           ENDIF
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      _EXCH_XYZ_R4(hFacC , myThid )

      DO bj=myByLo(myThid), myByHi(myThid)
       DO bi=myBxLo(myThid), myBxHi(myThid)
        DO K=1, Nz
         DO J=1,sNy
          DO I=1,sNx
           hFacW(I,J,K,bi,bj)=
     &       MIN(hFacC(I,J,K,bi,bj),hFacC(I-1,J,K,bi,bj))
           hFacS(I,J,K,bi,bj)=
     &       MIN(hFacC(I,J,K,bi,bj),hFacC(I,J-1,K,bi,bj))
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      _EXCH_XYZ_R4(hFacW , myThid )
      _EXCH_XYZ_R4(hFacS , myThid )

C--   Calculate recipricols of hFacC, hFacW and hFacS
      DO bj=myByLo(myThid), myByHi(myThid)
       DO bi=myBxLo(myThid), myBxHi(myThid)
        DO K=1, Nz
         DO J=1,sNy
          DO I=1,sNx
           rhFacC(I,J,K,bi,bj)=0. D0
           if (hFacC(I,J,K,bi,bj).ne.0.)
     &            rhFacC(I,J,K,bi,bj)=1. D0 /hFacC(I,J,K,bi,bj)
           rhFacW(I,J,K,bi,bj)=0. D0
           if (hFacW(I,J,K,bi,bj).ne.0.)
     &            rhFacW(I,J,K,bi,bj)=1. D0 /hFacW(I,J,K,bi,bj)
           rhFacS(I,J,K,bi,bj)=0. D0
           if (hFacS(I,J,K,bi,bj).ne.0.)
     &            rhFacS(I,J,K,bi,bj)=1. D0 /hFacS(I,J,K,bi,bj)
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      _EXCH_XYZ_R4(rhFacC , myThid )
      _EXCH_XYZ_R4(rhFacW , myThid )
      _EXCH_XYZ_R4(rhFacS , myThid )
C
      RETURN
      END
1	cnh	1.3	C $Header: /u/gcmpack/models/MITgcmUV/model/src/ini_depths.F,v 1.2 1998/04/24 02:05:41 cnh Exp $
2	cnh	1.1
3			#include "CPP_EEOPTIONS.h"
4
5			CStartOfInterface
6			SUBROUTINE INI_DEPTHS( myThid )
7			C /==========================================================\
8			C \| SUBROUTINE INI_DEPTHS \|
9			C \| o Initialise map of model depths \|
10			C \|==========================================================\|
11			C \| The depths of the bottom of the model is specified in \|
12			C \| terms of an XY map with one depth for each column of \|
13			C \| grid cells. Depths do not have to coincide with the \|
14			C \| model levels. The model's lopping algorithm makes it \|
15			C \| possible to represent arbitrary depths. \|
16			C \| The mode depths map also influences the models topology \|
17			C \| By default the model domain wraps around in X and Y. \|
18			C \| This default doubly periodic topology is "supressed" \|
19			C \| if a depth map is defined which closes off all wrap \|
20			C \| around flow. \|
21			C \==========================================================/
22
23			C === Global variables ===
24			#include "SIZE.h"
25			#include "EEPARAMS.h"
26			#include "PARAMS.h"
27			#include "GRID.h"
28
29			C == Routine arguments ==
30			C myThid - Number of this instance of INI_DEPTHS
31			INTEGER myThid
32			CEndOfInterface
33
34			C == Local variables ==
35			C xG, yG - Global coordinate location.
36			C zG
37			C zUpper - Work arrays for upper and lower
38			C zLower cell-face heights.
39			C phi - Temporary scalar
40			C iG, jG - Global coordinate index
41			C bi,bj - Loop counters
42			C zUpper - Temporary arrays holding z coordinates of
43			C zLower upper and lower faces.
44			C I,J,K
45			_RL xG, yG, zG
46			_RL phi
47			_RL zUpper(Nz), zLower(Nz)
48			INTEGER iG, jG
49			INTEGER bi, bj
50			INTEGER I, J, K
51
52			C For now set up a flat bottom box with doubly periodic topology.
53			C H is the basic variable from which other terms are derived. It
54			C is the term that would be set from an external file for a
55			C realistic problem.
56			_BARRIER
57			phi = 0. D0
58			DO K=1,Nz
59			phi = phi+delZ(K)
60			ENDDO
61			DO bj = myByLo(myThid), myByHi(myThid)
62			DO bi = myBxLo(myThid), myBxHi(myThid)
63			DO j=1,sNy
64			DO i=1,sNx
65			iG = myXGlobalLo-1+(bi-1)*sNx+I
66			jG = myYGlobalLo-1+(bj-1)*sNy+J
67			C Default depth of full domain
68			H(i,j,bi,bj) = phi
69			C Test for eastern edge
70			IF ( iG .EQ. nX ) H(i,j,bi,bj) = 0.
71			C Test for northern edge
72			IF ( jG .EQ. nY ) H(i,j,bi,bj) = 0.
73			C Island
74	cnh	1.3	IF ( iG .EQ. 1 .AND.
75			& jG .EQ. 24 ) H(i,j,bi,bj) = 0.75*phi
76	cnh	1.1	C IF ( iG .EQ. 1 .AND.
77			C & jG .EQ. 24 ) H(i,j,bi,bj) = 0.
78
79			ENDDO
80			ENDDO
81			ENDDO
82			ENDDO
83			DO bj = myByLo(myThid), myByHi(myThid)
84			DO bi = myBxLo(myThid), myBxHi(myThid)
85			DO J=1,sNy
86			DO I=1,sNx
87			C Inverse of depth
88			IF ( h(i,j,bi,bj) .EQ. 0. _d 0 ) THEN
89			rH(i,j,bi,bj) = 0. _d 0
90			ELSE
91			rH(i,j,bi,bj) = 1. _d 0 / h(i,j,bi,bj)
92			ENDIF
93			ENDDO
94			ENDDO
95			ENDDO
96			ENDDO
97			_EXCH_XY_R4( H, myThid )
98			_EXCH_XY_R4( rH, myThid )
99
100			C-- Now calculate "lopping" factors hFac.
101			zG = delZ(1)*0.5 D0
102			zFace(1) = 0
103			rDzC(1) = 2. _d 0/delZ(1)
104			DO K=1,Nz-1
105			saFac(K) = 1. D0
106			zC(K) = zG
107			zG = zG + (delZ(K)+delZ(K+1))*0.5 D0
108			zFace(K+1) = zFace(K)+delZ(K)
109			rDzC(K+1) = 2. _d 0/(delZ(K)+delZ(K+1))
110			ENDDO
111			zC(Nz) = zG
112			zFace(Nz+1) = zFace(Nz)+delZ(Nz)
113			DO K=1,Nz
114			zUpper(K) = zFace(K)
115			zLower(K) = zFace(K+1)
116			dzF(K) = delz(K)
117			rdzF(K) = 1.d0/dzF(K)
118			ENDDO
119			DO bj=myByLo(myThid), myByHi(myThid)
120			DO bi=myBxLo(myThid), myBxHi(myThid)
121			DO K=1, Nz
122			DO J=1,sNy
123			DO I=1,sNx
124			hFacC(I,J,K,bi,bj) = 1.
125			IF ( H(I,J,bi,bj) .LE. zUpper(K) ) THEN
126			C Below base of domain
127			hFacC(I,J,K,bi,bj) = 0.
128			ELSEIF ( H(I,J,bi,bj) .GT. zLower(K) ) THEN
129			C Base of domain is below this cell
130			hFacC(I,J,K,bi,bj) = 1.
131			ELSE
132			C Base of domain is in this cell
133			C Set hFac tp the fraction of the cell that is open.
134			phi = zUpper(K)-H(I,J,bi,bj)
135			hFacC(I,J,K,bi,bj) = phi/(zUpper(K)-zLower(K))
136			ENDIF
137			ENDDO
138			ENDDO
139			ENDDO
140			ENDDO
141			ENDDO
142			_EXCH_XYZ_R4(hFacC , myThid )
143
144			DO bj=myByLo(myThid), myByHi(myThid)
145			DO bi=myBxLo(myThid), myBxHi(myThid)
146			DO K=1, Nz
147			DO J=1,sNy
148			DO I=1,sNx
149			hFacW(I,J,K,bi,bj)=
150			& MIN(hFacC(I,J,K,bi,bj),hFacC(I-1,J,K,bi,bj))
151			hFacS(I,J,K,bi,bj)=
152			& MIN(hFacC(I,J,K,bi,bj),hFacC(I,J-1,K,bi,bj))
153			ENDDO
154			ENDDO
155			ENDDO
156			ENDDO
157			ENDDO
158			_EXCH_XYZ_R4(hFacW , myThid )
159			_EXCH_XYZ_R4(hFacS , myThid )
160
161			C-- Calculate recipricols of hFacC, hFacW and hFacS
162			DO bj=myByLo(myThid), myByHi(myThid)
163			DO bi=myBxLo(myThid), myBxHi(myThid)
164			DO K=1, Nz
165			DO J=1,sNy
166			DO I=1,sNx
167			rhFacC(I,J,K,bi,bj)=0. D0
168			if (hFacC(I,J,K,bi,bj).ne.0.)
169			& rhFacC(I,J,K,bi,bj)=1. D0 /hFacC(I,J,K,bi,bj)
170			rhFacW(I,J,K,bi,bj)=0. D0
171			if (hFacW(I,J,K,bi,bj).ne.0.)
172			& rhFacW(I,J,K,bi,bj)=1. D0 /hFacW(I,J,K,bi,bj)
173			rhFacS(I,J,K,bi,bj)=0. D0
174			if (hFacS(I,J,K,bi,bj).ne.0.)
175			& rhFacS(I,J,K,bi,bj)=1. D0 /hFacS(I,J,K,bi,bj)
176			ENDDO
177			ENDDO
178			ENDDO
179			ENDDO
180			ENDDO
181			_EXCH_XYZ_R4(rhFacC , myThid )
182			_EXCH_XYZ_R4(rhFacW , myThid )
183			_EXCH_XYZ_R4(rhFacS , myThid )
184			C
185			RETURN
186			END