/[MITgcm]/MITgcm/model/src/ini_masks_etc.F
ViewVC logotype

Contents of /MITgcm/model/src/ini_masks_etc.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.17 - (show annotations) (download)
Tue May 18 17:57:01 1999 UTC (25 years ago) by adcroft
Branch: MAIN
CVS Tags: checkpoint22, checkpoint23, checkpoint24, checkpoint25, checkpoint27, checkpoint26
Changes since 1.16: +3 -3 lines
Corrected calculation of depthinK(:,:,:). Instigated by Ralf.

1 C $Header: /u/gcmpack/models/MITgcmUV/model/src/ini_masks_etc.F,v 1.16 1999/05/05 18:32:34 adcroft Exp $
2
3 #include "CPP_OPTIONS.h"
4
5 CStartOfInterface
6 SUBROUTINE INI_MASKS_ETC( myThid )
7 C /==========================================================\
8 C | SUBROUTINE INI_MASKS_ETC |
9 C | o Initialise masks and topography factors |
10 C |==========================================================|
11 C | These arrays are used throughout the code and describe |
12 C | the topography of the domain through masks (0s and 1s) |
13 C | and fractional height factors (0<hFac<1). The latter |
14 C | distinguish between the lopped-cell and full-step |
15 C | topographic representations. |
16 C \==========================================================/
17 IMPLICIT NONE
18
19 C === Global variables ===
20 #include "SIZE.h"
21 #include "EEPARAMS.h"
22 #include "PARAMS.h"
23 #include "GRID.h"
24
25 C == Routine arguments ==
26 C myThid - Number of this instance of INI_MASKS_ETC
27 INTEGER myThid
28 CEndOfInterface
29
30 C == Local variables ==
31 C bi,bj - Loop counters
32 C I,J,K
33 INTEGER bi, bj
34 INTEGER I, J, K
35 #ifdef ALLOW_NONHYDROSTATIC
36 INTEGER Km1
37 _RL hFacUpper,hFacLower
38 #endif
39
40 C Calculate lopping factor hFacC
41 DO bj=myByLo(myThid), myByHi(myThid)
42 DO bi=myBxLo(myThid), myBxHi(myThid)
43 DO K=1, Nr
44 DO J=1,sNy
45 DO I=1,sNx
46 C Round depths within a small fraction of layer depth to that
47 C layer depth.
48 IF ( ABS(H(I,J,bi,bj)-rF(K)) .LT.
49 & 1. _d -6*ABS(rF(K)) .AND.
50 & ABS(H(I,J,bi,bj)-rF(K)) .LT.
51 & 1. _d -6*ABS(H(I,J,bi,bj)) )THEN
52 H(I,J,bi,bj) = rF(K)
53 ENDIF
54 IF ( H(I,J,bi,bj)*rkFac .GE. rF(K)*rkFac ) THEN
55 C Top of cell is below base of domain
56 hFacC(I,J,K,bi,bj) = 0.
57 ELSEIF ( H(I,J,bi,bj)*rkFac .LE. rF(K+1)*rkFac ) THEN
58 C Base of domain is below bottom of this cell
59 hFacC(I,J,K,bi,bj) = 1.
60 ELSE
61 C Base of domain is in this cell
62 C Set hFac to the fraction of the cell that is open.
63 hFacC(I,J,K,bi,bj) =
64 & (rF(K)*rkFac-H(I,J,bi,bj)*rkFac)*recip_drF(K)
65 ENDIF
66 C Impose minimum fraction
67 IF (hFacC(I,J,K,bi,bj).LT.hFacMin) THEN
68 IF (hFacC(I,J,K,bi,bj).LT.hFacMin*0.5) THEN
69 hFacC(I,J,K,bi,bj)=0.
70 ELSE
71 hFacC(I,J,K,bi,bj)=hFacMin
72 ENDIF
73 ENDIF
74 C Impose minimum size (dimensional)
75 IF (drF(k)*hFacC(I,J,K,bi,bj).LT.hFacMinDr) THEN
76 IF (drF(k)*hFacC(I,J,K,bi,bj).LT.hFacMinDr*0.5) THEN
77 hFacC(I,J,K,bi,bj)=0.
78 ELSE
79 hFacC(I,J,K,bi,bj)=hFacMinDr*recip_drF(k)
80 ENDIF
81 ENDIF
82 depthInK(i,j,bi,bj) = depthInK(i,j,bi,bj) + 1.
83 Crg & +hFacC(i,j,k,bi,bj)
84 ENDDO
85 ENDDO
86 ENDDO
87 ENDDO
88 ENDDO
89 _EXCH_XYZ_R4(hFacC , myThid )
90 _EXCH_XY_R4( depthInK, myThid )
91
92 CALL PLOT_FIELD_XYRS( depthInK,
93 & 'Model Depths K Index' , 1, myThid )
94
95 C Re-calculate depth of ocean, taking into account hFacC
96 DO bj=myByLo(myThid), myByHi(myThid)
97 DO bi=myBxLo(myThid), myBxHi(myThid)
98 DO J=1,sNy
99 DO I=1,sNx
100 H(I,J,bi,bj)=0.
101 DO K=1,Nr
102 H(I,J,bi,bj)=H(I,J,bi,bj)-
103 & rkFac*drF(k)*hFacC(I,J,K,bi,bj)
104 ENDDO
105 ENDDO
106 ENDDO
107 ENDDO
108 ENDDO
109 _EXCH_XY_R4(H , myThid )
110 CALL WRITE_FLD_XY_RS( 'Depth',' ',H,0,myThid)
111 C CALL MDSWRITEFIELD( 'Depth', writeBinaryPrec, .TRUE.,
112 C & 'RS', 1, H, 1, -1, myThid )
113
114 C Calculate quantities derived from XY depth map
115 DO bj = myByLo(myThid), myByHi(myThid)
116 DO bi = myBxLo(myThid), myBxHi(myThid)
117 DO J=1,sNy
118 DO I=1,sNx
119 C Inverse of depth
120 IF ( h(i,j,bi,bj) .EQ. 0. _d 0 ) THEN
121 recip_H(i,j,bi,bj) = 0. _d 0
122 ELSE
123 recip_H(i,j,bi,bj) = 1. _d 0 / abs( H(i,j,bi,bj) )
124 ENDIF
125 depthInK(i,j,bi,bj) = 0.
126 ENDDO
127 ENDDO
128 ENDDO
129 ENDDO
130 _EXCH_XY_R4( recip_H, myThid )
131
132 C hFacW and hFacS (at U and V points)
133 DO bj=myByLo(myThid), myByHi(myThid)
134 DO bi=myBxLo(myThid), myBxHi(myThid)
135 DO K=1, Nr
136 DO J=1,sNy
137 DO I=1,sNx
138 hFacW(I,J,K,bi,bj)=
139 & MIN(hFacC(I,J,K,bi,bj),hFacC(I-1,J,K,bi,bj))
140 hFacS(I,J,K,bi,bj)=
141 & MIN(hFacC(I,J,K,bi,bj),hFacC(I,J-1,K,bi,bj))
142 ENDDO
143 ENDDO
144 ENDDO
145 ENDDO
146 ENDDO
147 _EXCH_XYZ_R4(hFacW , myThid )
148 _EXCH_XYZ_R4(hFacS , myThid )
149
150 C Masks and reciprocals of hFac[CWS]
151 DO bj = myByLo(myThid), myByHi(myThid)
152 DO bi = myBxLo(myThid), myBxHi(myThid)
153 DO K=1,Nr
154 DO J=1,sNy
155 DO I=1,sNx
156 IF (HFacC(I,J,K,bi,bj) .NE. 0. _d 0 ) THEN
157 recip_HFacC(I,J,K,bi,bj) = 1. _d 0 / HFacC(I,J,K,bi,bj)
158 ELSE
159 recip_HFacC(I,J,K,bi,bj) = 0. _d 0
160 ENDIF
161 IF (HFacW(I,J,K,bi,bj) .NE. 0. _d 0 ) THEN
162 recip_HFacW(I,J,K,bi,bj) = 1. _d 0 / HFacW(I,J,K,bi,bj)
163 maskW(I,J,K,bi,bj) = 1. _d 0
164 ELSE
165 recip_HFacW(I,J,K,bi,bj) = 0. _d 0
166 maskW(I,J,K,bi,bj) = 0.0 _d 0
167 ENDIF
168 IF (HFacS(I,J,K,bi,bj) .NE. 0. _d 0 ) THEN
169 recip_HFacS(I,J,K,bi,bj) = 1. _d 0 / HFacS(I,J,K,bi,bj)
170 maskS(I,J,K,bi,bj) = 1. _d 0
171 ELSE
172 recip_HFacS(I,J,K,bi,bj) = 0. _d 0
173 maskS(I,J,K,bi,bj) = 0. _d 0
174 ENDIF
175 ENDDO
176 ENDDO
177 ENDDO
178 ENDDO
179 ENDDO
180 _EXCH_XYZ_R4(recip_HFacC , myThid )
181 _EXCH_XYZ_R4(recip_HFacW , myThid )
182 _EXCH_XYZ_R4(recip_HFacS , myThid )
183 _EXCH_XYZ_R4(maskW , myThid )
184 _EXCH_XYZ_R4(maskS , myThid )
185
186 C Calculate recipricols grid lengths
187 DO bj = myByLo(myThid), myByHi(myThid)
188 DO bi = myBxLo(myThid), myBxHi(myThid)
189 DO J=1,sNy
190 DO I=1,sNx
191 recip_dxG(I,J,bi,bj)=1.d0/dxG(I,J,bi,bj)
192 recip_dyG(I,J,bi,bj)=1.d0/dyG(I,J,bi,bj)
193 recip_dxC(I,J,bi,bj)=1.d0/dxC(I,J,bi,bj)
194 recip_dyC(I,J,bi,bj)=1.d0/dyC(I,J,bi,bj)
195 recip_dxF(I,J,bi,bj)=1.d0/dxF(I,J,bi,bj)
196 recip_dyF(I,J,bi,bj)=1.d0/dyF(I,J,bi,bj)
197 recip_dxV(I,J,bi,bj)=1.d0/dxV(I,J,bi,bj)
198 recip_dyU(I,J,bi,bj)=1.d0/dyU(I,J,bi,bj)
199 ENDDO
200 ENDDO
201 ENDDO
202 ENDDO
203 _EXCH_XY_R4(recip_dxG, myThid )
204 _EXCH_XY_R4(recip_dyG, myThid )
205 _EXCH_XY_R4(recip_dxC, myThid )
206 _EXCH_XY_R4(recip_dyC, myThid )
207 _EXCH_XY_R4(recip_dxF, myThid )
208 _EXCH_XY_R4(recip_dyF, myThid )
209 _EXCH_XY_R4(recip_dxV, myThid )
210 _EXCH_XY_R4(recip_dyU, myThid )
211
212 #ifdef ALLOW_NONHYDROSTATIC
213 C-- Calculate the reciprocal hfac distance/volume for W cells
214 DO bj = myByLo(myThid), myByHi(myThid)
215 DO bi = myBxLo(myThid), myBxHi(myThid)
216 DO K=1,Nr
217 Km1=max(K-1,1)
218 hFacUpper=drF(Km1)/(drF(Km1)+drF(K))
219 IF (Km1.EQ.K) hFacUpper=0.
220 hFacLower=drF(K)/(drF(Km1)+drF(K))
221 DO J=1,sNy
222 DO I=1,sNx
223 IF (hFacC(I,J,K,bi,bj).NE.0.) THEN
224 IF (hFacC(I,J,K,bi,bj).LE.0.5) THEN
225 recip_hFacU(I,J,K,bi,bj)=
226 & hFacUpper+hFacLower*hFacC(I,J,K,bi,bj)
227 ELSE
228 recip_hFacU(I,J,K,bi,bj)=1.
229 ENDIF
230 ELSE
231 recip_hFacU(I,J,K,bi,bj)=0.
232 ENDIF
233 IF (recip_hFacU(I,J,K,bi,bj).NE.0.)
234 & recip_hFacU(I,J,K,bi,bj)=1./recip_hFacU(I,J,K,bi,bj)
235 ENDDO
236 ENDDO
237 ENDDO
238 ENDDO
239 ENDDO
240 _EXCH_XY_R4(recip_hFacU, myThid )
241 #endif
242 C
243 RETURN
244 END

  ViewVC Help
Powered by ViewVC 1.1.22