/[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.18 - (show annotations) (download)
Tue Apr 11 13:39:09 2000 UTC (24 years, 1 month ago) by adcroft
Branch: MAIN
CVS Tags: checkpoint28, checkpoint29, branch-atmos-merge-start, checkpoint33, checkpoint32, checkpoint31, checkpoint30, checkpoint34, branch-atmos-merge-phase5, branch-atmos-merge-phase4, branch-atmos-merge-phase6, branch-atmos-merge-phase1, branch-atmos-merge-phase3, branch-atmos-merge-phase2
Branch point for: branch-atmos-merge
Changes since 1.17: +7 -13 lines 
Modified the "lopping" procedure for determining the fraction cell
height.

The original formulation first filter the model depth using
the hFacMin parameter. Any fractional cells less than hFacMin
were rounded to the nearer of 0 or hFacMin. Then the depths were
filtered using hFacMinDr where the actual size was tested.
This meant that some cells could be rounded twice and also
that if hFacMinDr was specified as larger than the nominal
thickness that the fractional cell thickness could end up
larger than 1. Clearly a bad idea (?) but to cause this one
would have to specify the hFacMin and hFacMinDr parameters in
a inconsistent way.

The new method combines hFacMin and hFacMinDr first and then
does a single filtering of model depths. If hFacMinDr is
greater than the nominal layer thickness, the later is used
and the fractional cell thicknesses are set to 1. Otherwise,
layers thicknesses are rounded to using the larger of
hFacMin and hFacMinDr/drF(K).

Using the new method, setting hFacMin=0 and hFacMin=100
means that lopping only occurs for cell thicker than 100
and the minimum cell thickness of lopped cells will be 100.
This, I think, is a much cleaner specification. AJA.
1 C $Header: /u/gcmpack/models/MITgcmUV/model/src/ini_masks_etc.F,v 1.17 1999/05/18 17:57:01 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 _RL hFacTmp
36 #ifdef ALLOW_NONHYDROSTATIC
37 INTEGER Km1
38 _RL hFacUpper,hFacLower
39 #endif
40
41 C Calculate lopping factor hFacC
42 DO bj=myByLo(myThid), myByHi(myThid)
43 DO bi=myBxLo(myThid), myBxHi(myThid)
44 DO K=1, Nr
45 DO J=1,sNy
46 DO I=1,sNx
47 C Round depths within a small fraction of layer depth to that
48 C layer depth.
49 IF ( ABS(H(I,J,bi,bj)-rF(K)) .LT.
50 & 1. _d -6*ABS(rF(K)) .AND.
51 & ABS(H(I,J,bi,bj)-rF(K)) .LT.
52 & 1. _d -6*ABS(H(I,J,bi,bj)) )THEN
53 H(I,J,bi,bj) = rF(K)
54 ENDIF
55 IF ( H(I,J,bi,bj)*rkFac .GE. rF(K)*rkFac ) THEN
56 C Top of cell is below base of domain
57 hFacC(I,J,K,bi,bj) = 0.
58 ELSEIF ( H(I,J,bi,bj)*rkFac .LE. rF(K+1)*rkFac ) THEN
59 C Base of domain is below bottom of this cell
60 hFacC(I,J,K,bi,bj) = 1.
61 ELSE
62 C Base of domain is in this cell
63 C Set hFac to the fraction of the cell that is open.
64 hFacC(I,J,K,bi,bj) =
65 & (rF(K)*rkFac-H(I,J,bi,bj)*rkFac)*recip_drF(K)
66 ENDIF
67 C Impose minimum fraction and/or size
68 hFacTmp=max( hFacMin , min(hFacMinDr*recip_drF(k),1.) )
69 IF (hFacC(I,J,K,bi,bj).LT.hFacTmp) THEN
70 IF (hFacC(I,J,K,bi,bj).LT.hFacTmp*0.5) THEN
71 hFacC(I,J,K,bi,bj)=0.
72 ELSE
73 hFacC(I,J,K,bi,bj)=hFacTmp
74 ENDIF
75 ENDIF
76 depthInK(i,j,bi,bj) = depthInK(i,j,bi,bj) + 1.
77 Crg & +hFacC(i,j,k,bi,bj)
78 ENDDO
79 ENDDO
80 ENDDO
81 ENDDO
82 ENDDO
83 _EXCH_XYZ_R4(hFacC , myThid )
84 _EXCH_XY_R4( depthInK, myThid )
85
86 CALL PLOT_FIELD_XYRS( depthInK,
87 & 'Model Depths K Index' , 1, myThid )
88
89 C Re-calculate depth of ocean, taking into account hFacC
90 DO bj=myByLo(myThid), myByHi(myThid)
91 DO bi=myBxLo(myThid), myBxHi(myThid)
92 DO J=1,sNy
93 DO I=1,sNx
94 H(I,J,bi,bj)=0.
95 DO K=1,Nr
96 H(I,J,bi,bj)=H(I,J,bi,bj)-
97 & rkFac*drF(k)*hFacC(I,J,K,bi,bj)
98 ENDDO
99 ENDDO
100 ENDDO
101 ENDDO
102 ENDDO
103 _EXCH_XY_R4(H , myThid )
104 CALL WRITE_FLD_XY_RS( 'Depth',' ',H,0,myThid)
105 C CALL MDSWRITEFIELD( 'Depth', writeBinaryPrec, .TRUE.,
106 C & 'RS', 1, H, 1, -1, myThid )
107
108 C Calculate quantities derived from XY depth map
109 DO bj = myByLo(myThid), myByHi(myThid)
110 DO bi = myBxLo(myThid), myBxHi(myThid)
111 DO J=1,sNy
112 DO I=1,sNx
113 C Inverse of depth
114 IF ( h(i,j,bi,bj) .EQ. 0. _d 0 ) THEN
115 recip_H(i,j,bi,bj) = 0. _d 0
116 ELSE
117 recip_H(i,j,bi,bj) = 1. _d 0 / abs( H(i,j,bi,bj) )
118 ENDIF
119 depthInK(i,j,bi,bj) = 0.
120 ENDDO
121 ENDDO
122 ENDDO
123 ENDDO
124 _EXCH_XY_R4( recip_H, myThid )
125
126 C hFacW and hFacS (at U and V points)
127 DO bj=myByLo(myThid), myByHi(myThid)
128 DO bi=myBxLo(myThid), myBxHi(myThid)
129 DO K=1, Nr
130 DO J=1,sNy
131 DO I=1,sNx
132 hFacW(I,J,K,bi,bj)=
133 & MIN(hFacC(I,J,K,bi,bj),hFacC(I-1,J,K,bi,bj))
134 hFacS(I,J,K,bi,bj)=
135 & MIN(hFacC(I,J,K,bi,bj),hFacC(I,J-1,K,bi,bj))
136 ENDDO
137 ENDDO
138 ENDDO
139 ENDDO
140 ENDDO
141 _EXCH_XYZ_R4(hFacW , myThid )
142 _EXCH_XYZ_R4(hFacS , myThid )
143
144 C Masks and reciprocals of hFac[CWS]
145 DO bj = myByLo(myThid), myByHi(myThid)
146 DO bi = myBxLo(myThid), myBxHi(myThid)
147 DO K=1,Nr
148 DO J=1,sNy
149 DO I=1,sNx
150 IF (HFacC(I,J,K,bi,bj) .NE. 0. _d 0 ) THEN
151 recip_HFacC(I,J,K,bi,bj) = 1. _d 0 / HFacC(I,J,K,bi,bj)
152 ELSE
153 recip_HFacC(I,J,K,bi,bj) = 0. _d 0
154 ENDIF
155 IF (HFacW(I,J,K,bi,bj) .NE. 0. _d 0 ) THEN
156 recip_HFacW(I,J,K,bi,bj) = 1. _d 0 / HFacW(I,J,K,bi,bj)
157 maskW(I,J,K,bi,bj) = 1. _d 0
158 ELSE
159 recip_HFacW(I,J,K,bi,bj) = 0. _d 0
160 maskW(I,J,K,bi,bj) = 0.0 _d 0
161 ENDIF
162 IF (HFacS(I,J,K,bi,bj) .NE. 0. _d 0 ) THEN
163 recip_HFacS(I,J,K,bi,bj) = 1. _d 0 / HFacS(I,J,K,bi,bj)
164 maskS(I,J,K,bi,bj) = 1. _d 0
165 ELSE
166 recip_HFacS(I,J,K,bi,bj) = 0. _d 0
167 maskS(I,J,K,bi,bj) = 0. _d 0
168 ENDIF
169 ENDDO
170 ENDDO
171 ENDDO
172 ENDDO
173 ENDDO
174 _EXCH_XYZ_R4(recip_HFacC , myThid )
175 _EXCH_XYZ_R4(recip_HFacW , myThid )
176 _EXCH_XYZ_R4(recip_HFacS , myThid )
177 _EXCH_XYZ_R4(maskW , myThid )
178 _EXCH_XYZ_R4(maskS , myThid )
179
180 C Calculate recipricols grid lengths
181 DO bj = myByLo(myThid), myByHi(myThid)
182 DO bi = myBxLo(myThid), myBxHi(myThid)
183 DO J=1,sNy
184 DO I=1,sNx
185 recip_dxG(I,J,bi,bj)=1.d0/dxG(I,J,bi,bj)
186 recip_dyG(I,J,bi,bj)=1.d0/dyG(I,J,bi,bj)
187 recip_dxC(I,J,bi,bj)=1.d0/dxC(I,J,bi,bj)
188 recip_dyC(I,J,bi,bj)=1.d0/dyC(I,J,bi,bj)
189 recip_dxF(I,J,bi,bj)=1.d0/dxF(I,J,bi,bj)
190 recip_dyF(I,J,bi,bj)=1.d0/dyF(I,J,bi,bj)
191 recip_dxV(I,J,bi,bj)=1.d0/dxV(I,J,bi,bj)
192 recip_dyU(I,J,bi,bj)=1.d0/dyU(I,J,bi,bj)
193 ENDDO
194 ENDDO
195 ENDDO
196 ENDDO
197 _EXCH_XY_R4(recip_dxG, myThid )
198 _EXCH_XY_R4(recip_dyG, myThid )
199 _EXCH_XY_R4(recip_dxC, myThid )
200 _EXCH_XY_R4(recip_dyC, myThid )
201 _EXCH_XY_R4(recip_dxF, myThid )
202 _EXCH_XY_R4(recip_dyF, myThid )
203 _EXCH_XY_R4(recip_dxV, myThid )
204 _EXCH_XY_R4(recip_dyU, myThid )
205
206 #ifdef ALLOW_NONHYDROSTATIC
207 C-- Calculate the reciprocal hfac distance/volume for W cells
208 DO bj = myByLo(myThid), myByHi(myThid)
209 DO bi = myBxLo(myThid), myBxHi(myThid)
210 DO K=1,Nr
211 Km1=max(K-1,1)
212 hFacUpper=drF(Km1)/(drF(Km1)+drF(K))
213 IF (Km1.EQ.K) hFacUpper=0.
214 hFacLower=drF(K)/(drF(Km1)+drF(K))
215 DO J=1,sNy
216 DO I=1,sNx
217 IF (hFacC(I,J,K,bi,bj).NE.0.) THEN
218 IF (hFacC(I,J,K,bi,bj).LE.0.5) THEN
219 recip_hFacU(I,J,K,bi,bj)=
220 & hFacUpper+hFacLower*hFacC(I,J,K,bi,bj)
221 ELSE
222 recip_hFacU(I,J,K,bi,bj)=1.
223 ENDIF
224 ELSE
225 recip_hFacU(I,J,K,bi,bj)=0.
226 ENDIF
227 IF (recip_hFacU(I,J,K,bi,bj).NE.0.)
228 & recip_hFacU(I,J,K,bi,bj)=1./recip_hFacU(I,J,K,bi,bj)
229 ENDDO
230 ENDDO
231 ENDDO
232 ENDDO
233 ENDDO
234 _EXCH_XY_R4(recip_hFacU, myThid )
235 #endif
236 C
237 RETURN
238 END
  ViewVC Help
Powered by ViewVC 1.1.22