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

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

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

Revision 1.18 - (hide 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 adcroft 1.18 C $Header: /u/gcmpack/models/MITgcmUV/model/src/ini_masks_etc.F,v 1.17 1999/05/18 17:57:01 adcroft Exp $
2 adcroft 1.1
3 cnh 1.11 #include "CPP_OPTIONS.h"
4 adcroft 1.1
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 adcroft 1.13 IMPLICIT NONE
18 adcroft 1.1
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 cnh 1.6 C myThid - Number of this instance of INI_MASKS_ETC
27 adcroft 1.1 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 adcroft 1.18 _RL hFacTmp
36 adcroft 1.15 #ifdef ALLOW_NONHYDROSTATIC
37     INTEGER Km1
38     _RL hFacUpper,hFacLower
39     #endif
40 adcroft 1.1
41 adcroft 1.2 C Calculate lopping factor hFacC
42     DO bj=myByLo(myThid), myByHi(myThid)
43     DO bi=myBxLo(myThid), myBxHi(myThid)
44 cnh 1.4 DO K=1, Nr
45 adcroft 1.2 DO J=1,sNy
46     DO I=1,sNx
47 cnh 1.7 C Round depths within a small fraction of layer depth to that
48     C layer depth.
49 cnh 1.9 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 cnh 1.7 H(I,J,bi,bj) = rF(K)
54     ENDIF
55 cnh 1.6 IF ( H(I,J,bi,bj)*rkFac .GE. rF(K)*rkFac ) THEN
56 adcroft 1.3 C Top of cell is below base of domain
57 adcroft 1.2 hFacC(I,J,K,bi,bj) = 0.
58 cnh 1.6 ELSEIF ( H(I,J,bi,bj)*rkFac .LE. rF(K+1)*rkFac ) THEN
59 adcroft 1.3 C Base of domain is below bottom of this cell
60 adcroft 1.2 hFacC(I,J,K,bi,bj) = 1.
61     ELSE
62     C Base of domain is in this cell
63 adcroft 1.3 C Set hFac to the fraction of the cell that is open.
64 adcroft 1.15 hFacC(I,J,K,bi,bj) =
65     & (rF(K)*rkFac-H(I,J,bi,bj)*rkFac)*recip_drF(K)
66 adcroft 1.3 ENDIF
67 adcroft 1.18 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 adcroft 1.3 hFacC(I,J,K,bi,bj)=0.
72     ELSE
73 adcroft 1.18 hFacC(I,J,K,bi,bj)=hFacTmp
74 adcroft 1.3 ENDIF
75 adcroft 1.2 ENDIF
76 adcroft 1.17 depthInK(i,j,bi,bj) = depthInK(i,j,bi,bj) + 1.
77     Crg & +hFacC(i,j,k,bi,bj)
78 adcroft 1.2 ENDDO
79     ENDDO
80     ENDDO
81     ENDDO
82     ENDDO
83     _EXCH_XYZ_R4(hFacC , myThid )
84 cnh 1.7 _EXCH_XY_R4( depthInK, myThid )
85    
86 cnh 1.9 CALL PLOT_FIELD_XYRS( depthInK,
87     & 'Model Depths K Index' , 1, myThid )
88 adcroft 1.16
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 adcroft 1.1
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 cnh 1.4 DO K=1, Nr
130 adcroft 1.1 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 cnh 1.4 DO K=1,Nr
148 adcroft 1.1 DO J=1,sNy
149     DO I=1,sNx
150 cnh 1.10 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 adcroft 1.1 ELSE
153 cnh 1.10 recip_HFacC(I,J,K,bi,bj) = 0. _d 0
154 adcroft 1.1 ENDIF
155 cnh 1.10 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 adcroft 1.1 ELSE
159 cnh 1.10 recip_HFacW(I,J,K,bi,bj) = 0. _d 0
160     maskW(I,J,K,bi,bj) = 0.0 _d 0
161 adcroft 1.1 ENDIF
162 cnh 1.10 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 adcroft 1.1 ELSE
166 cnh 1.10 recip_HFacS(I,J,K,bi,bj) = 0. _d 0
167     maskS(I,J,K,bi,bj) = 0. _d 0
168 adcroft 1.1 ENDIF
169     ENDDO
170     ENDDO
171     ENDDO
172     ENDDO
173     ENDDO
174 cnh 1.4 _EXCH_XYZ_R4(recip_HFacC , myThid )
175     _EXCH_XYZ_R4(recip_HFacW , myThid )
176     _EXCH_XYZ_R4(recip_HFacS , myThid )
177 adcroft 1.1 _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 cnh 1.4 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 adcroft 1.1 ENDDO
194     ENDDO
195     ENDDO
196     ENDDO
197 cnh 1.4 _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 adcroft 1.1
206 adcroft 1.15 #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 adcroft 1.1 C
237     RETURN
238     END
  ViewVC Help
Powered by ViewVC 1.1.22