/[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.20 - (hide annotations) (download)
Sun Feb 4 14:38:47 2001 UTC (23 years, 3 months ago) by cnh
Branch: MAIN
CVS Tags: checkpoint38, c37_adj, checkpoint39, checkpoint37, checkpoint36, checkpoint35
Branch point for: pre38
Changes since 1.19: +2 -1 lines
Made sure each .F and .h file had
the CVS keywords Header and Name at its start.
Most had header but very few currently have Name, so
lots of changes!

1 cnh 1.20 C $Header: /u/gcmpack/models/MITgcmUV/model/src/ini_masks_etc.F,v 1.19 2001/02/02 21:04:48 adcroft Exp $
2     C $Name: $
3 adcroft 1.1
4 cnh 1.11 #include "CPP_OPTIONS.h"
5 adcroft 1.1
6     CStartOfInterface
7     SUBROUTINE INI_MASKS_ETC( myThid )
8     C /==========================================================\
9     C | SUBROUTINE INI_MASKS_ETC |
10     C | o Initialise masks and topography factors |
11     C |==========================================================|
12     C | These arrays are used throughout the code and describe |
13     C | the topography of the domain through masks (0s and 1s) |
14     C | and fractional height factors (0<hFac<1). The latter |
15     C | distinguish between the lopped-cell and full-step |
16     C | topographic representations. |
17     C \==========================================================/
18 adcroft 1.13 IMPLICIT NONE
19 adcroft 1.1
20     C === Global variables ===
21     #include "SIZE.h"
22     #include "EEPARAMS.h"
23     #include "PARAMS.h"
24     #include "GRID.h"
25    
26     C == Routine arguments ==
27 cnh 1.6 C myThid - Number of this instance of INI_MASKS_ETC
28 adcroft 1.1 INTEGER myThid
29     CEndOfInterface
30    
31     C == Local variables ==
32     C bi,bj - Loop counters
33     C I,J,K
34     INTEGER bi, bj
35     INTEGER I, J, K
36 adcroft 1.19 INTEGER kadj_rf, klev_noH
37 adcroft 1.15 #ifdef ALLOW_NONHYDROSTATIC
38     INTEGER Km1
39     _RL hFacUpper,hFacLower
40     #endif
41 adcroft 1.19 _RL hFacCtmp,topo_rkfac
42     _RL hFacMnSz
43    
44     IF (groundAtK1) THEN
45     topo_rkfac = -rkFac
46     kadj_rf = 1
47     klev_noH = Nr+1
48     ELSE
49     topo_rkfac = rkFac
50     kadj_rf = 0
51     klev_noH = 1
52     ENDIF
53 adcroft 1.1
54 adcroft 1.2 C Calculate lopping factor hFacC
55     DO bj=myByLo(myThid), myByHi(myThid)
56     DO bi=myBxLo(myThid), myBxHi(myThid)
57 cnh 1.4 DO K=1, Nr
58 adcroft 1.19 DO J=1-Oly,sNy+Oly
59     DO I=1-Olx,sNx+Olx
60     c IF (groundAtK1) THEN
61     C o Non-dimensional distance between grid boundary and model depth
62     C for case with "ground" at K=1 (i.e. like a good atmos model)
63     C e.g. rkfac=+1 => dR/dk<0 (eg. R=p): hFacCtmp=(H(x,y)-rF(k))/drF(K)
64     C e.g. rkfac=-1 => dR/dk>0 (eg. R=z): hFacCtmp=(rF(K)-H(x,y))/drF(K)
65     c hFacCtmp=rkFac*(H(I,J,bi,bj)-rF(K+1))*recip_drF(K)
66     c ELSE
67     C o Non-dimensional distance between grid boundary and model depth
68     C for case with "ground" at K=Nr (i.e. like original ocean model)
69     C e.g. rkfac=+1 => dR/dk<0 (eg. R=z): hFacCtmp=(rF(K)-H(x,y))/drF(K)
70     C e.g. rkfac=-1 => dR/dk>0 (eg. R=p): hFacCtmp=(H(x,y)-rF(k))/drF(K)
71     c hFacCtmp=rkFac*(rF(K)-H(I,J,bi,bj))*recip_drF(K)
72     c ENDIF
73     hFacCtmp=topo_rkfac*(rF(K+kadj_rf)-H(I,J,bi,bj))*recip_drF(K)
74     C o Select between, closed, open or partial (0,1,0-1)
75     hFacCtmp=min( max( hFacCtmp, 0. _d 0) , 1. _d 0)
76     C o And there we have it, the fractional open cell volume
77     hFacC(I,J,K,bi,bj)=hFacCtmp
78     C o Impose minimum fraction and/or size (dimensional)
79     hFacMnSz=max( hFacMin , min(hFacMinDr*recip_drF(k),1. _d 0) )
80     IF (hFacC(I,J,K,bi,bj).LT.hFacMnSz) THEN
81     IF (hFacC(I,J,K,bi,bj).LT.hFacMnSz*0.5) THEN
82 adcroft 1.3 hFacC(I,J,K,bi,bj)=0.
83     ELSE
84 adcroft 1.19 hFacC(I,J,K,bi,bj)=hFacMnSz
85 adcroft 1.3 ENDIF
86 adcroft 1.2 ENDIF
87 adcroft 1.19 IF (hFacC(I,J,K,bi,bj).NE.0.)
88     & depthInK(i,j,bi,bj) = depthInK(i,j,bi,bj) + 1.
89 adcroft 1.2 ENDDO
90     ENDDO
91     ENDDO
92     ENDDO
93     ENDDO
94 adcroft 1.19 C _EXCH_XYZ_R4(hFacC , myThid )
95     C _EXCH_XY_R4( depthInK, myThid )
96 cnh 1.7
97 cnh 1.9 CALL PLOT_FIELD_XYRS( depthInK,
98     & 'Model Depths K Index' , 1, myThid )
99 adcroft 1.16
100     C Re-calculate depth of ocean, taking into account hFacC
101     DO bj=myByLo(myThid), myByHi(myThid)
102     DO bi=myBxLo(myThid), myBxHi(myThid)
103 adcroft 1.19 DO J=1-Oly,sNy+Oly
104     DO I=1-Olx,sNx+Olx
105     H(I,J,bi,bj)=rF(klev_noH)
106 adcroft 1.16 DO K=1,Nr
107     H(I,J,bi,bj)=H(I,J,bi,bj)-
108 adcroft 1.19 & topo_rkFac*drF(k)*hFacC(I,J,K,bi,bj)
109 adcroft 1.16 ENDDO
110     ENDDO
111     ENDDO
112     ENDDO
113     ENDDO
114 adcroft 1.19 C _EXCH_XY_R4(H , myThid )
115     CALL PLOT_FIELD_XYRS(H,'Model depths (ini_masks_etc)',1,myThid)
116 adcroft 1.16 CALL WRITE_FLD_XY_RS( 'Depth',' ',H,0,myThid)
117 adcroft 1.19 CALL WRITE_FLD_XYZ_RS( 'hFacC',' ',hFacC,0,myThid)
118 adcroft 1.16 C CALL MDSWRITEFIELD( 'Depth', writeBinaryPrec, .TRUE.,
119     C & 'RS', 1, H, 1, -1, myThid )
120    
121     C Calculate quantities derived from XY depth map
122     DO bj = myByLo(myThid), myByHi(myThid)
123     DO bi = myBxLo(myThid), myBxHi(myThid)
124 adcroft 1.19 DO J=1-Oly,sNy+Oly
125     DO I=1-Olx,sNx+Olx
126 adcroft 1.16 C Inverse of depth
127 adcroft 1.19 IF ( h(i,j,bi,bj) .EQ. 0. ) THEN
128     recip_H(i,j,bi,bj) = 0.
129 adcroft 1.16 ELSE
130 adcroft 1.19 recip_H(i,j,bi,bj) = 1. / abs( H(i,j,bi,bj) )
131 adcroft 1.16 ENDIF
132     depthInK(i,j,bi,bj) = 0.
133     ENDDO
134     ENDDO
135     ENDDO
136     ENDDO
137 adcroft 1.19 C _EXCH_XY_R4( recip_H, myThid )
138 adcroft 1.1
139     C hFacW and hFacS (at U and V points)
140     DO bj=myByLo(myThid), myByHi(myThid)
141     DO bi=myBxLo(myThid), myBxHi(myThid)
142 cnh 1.4 DO K=1, Nr
143 adcroft 1.1 DO J=1,sNy
144     DO I=1,sNx
145     hFacW(I,J,K,bi,bj)=
146     & MIN(hFacC(I,J,K,bi,bj),hFacC(I-1,J,K,bi,bj))
147     hFacS(I,J,K,bi,bj)=
148     & MIN(hFacC(I,J,K,bi,bj),hFacC(I,J-1,K,bi,bj))
149     ENDDO
150     ENDDO
151     ENDDO
152     ENDDO
153     ENDDO
154     _EXCH_XYZ_R4(hFacW , myThid )
155     _EXCH_XYZ_R4(hFacS , myThid )
156 adcroft 1.19 C Re-do hFacW and hFacS (at U and V points)
157     DO bj=myByLo(myThid), myByHi(myThid)
158     DO bi=myBxLo(myThid), myBxHi(myThid)
159     DO K=1, Nr
160     DO J=1-Oly,sNy+Oly
161     DO I=1-Olx+1,sNx+Olx ! Note range
162     hFacW(I,J,K,bi,bj)=
163     & MIN(hFacC(I,J,K,bi,bj),hFacC(I-1,J,K,bi,bj))
164     ENDDO
165     ENDDO
166     DO I=1-Olx,sNx+Olx
167     DO J=1-Oly+1,sNy+Oly ! Note range
168     hFacS(I,J,K,bi,bj)=
169     & MIN(hFacC(I,J,K,bi,bj),hFacC(I,J-1,K,bi,bj))
170     ENDDO
171     ENDDO
172     ENDDO
173     ENDDO
174     ENDDO
175    
176     CALL PLOT_FIELD_XYZRS( hFacC, 'hFacC' , Nr, 1, myThid )
177     CALL PLOT_FIELD_XYZRS( hFacW, 'hFacW' , Nr, 1, myThid )
178     CALL PLOT_FIELD_XYZRS( hFacS, 'hFacS' , Nr, 1, myThid )
179 adcroft 1.1
180     C Masks and reciprocals of hFac[CWS]
181     DO bj = myByLo(myThid), myByHi(myThid)
182     DO bi = myBxLo(myThid), myBxHi(myThid)
183 cnh 1.4 DO K=1,Nr
184 adcroft 1.19 DO J=1-Oly,sNy+Oly
185     DO I=1-Olx,sNx+Olx
186     IF (HFacC(I,J,K,bi,bj) .NE. 0. ) THEN
187     recip_HFacC(I,J,K,bi,bj) = 1. / HFacC(I,J,K,bi,bj)
188 adcroft 1.1 ELSE
189 adcroft 1.19 recip_HFacC(I,J,K,bi,bj) = 0.
190 adcroft 1.1 ENDIF
191 adcroft 1.19 IF (HFacW(I,J,K,bi,bj) .NE. 0. ) THEN
192     recip_HFacW(I,J,K,bi,bj) = 1. / HFacW(I,J,K,bi,bj)
193     maskW(I,J,K,bi,bj) = 1.
194 adcroft 1.1 ELSE
195 adcroft 1.19 recip_HFacW(I,J,K,bi,bj) = 0.
196     maskW(I,J,K,bi,bj) = 0.
197 adcroft 1.1 ENDIF
198 adcroft 1.19 IF (HFacS(I,J,K,bi,bj) .NE. 0. ) THEN
199     recip_HFacS(I,J,K,bi,bj) = 1. / HFacS(I,J,K,bi,bj)
200     maskS(I,J,K,bi,bj) = 1.
201 adcroft 1.1 ELSE
202 adcroft 1.19 recip_HFacS(I,J,K,bi,bj) = 0.
203     maskS(I,J,K,bi,bj) = 0.
204 adcroft 1.1 ENDIF
205     ENDDO
206     ENDDO
207     ENDDO
208     ENDDO
209     ENDDO
210 adcroft 1.19 C _EXCH_XYZ_R4(recip_HFacC , myThid )
211     C _EXCH_XYZ_R4(recip_HFacW , myThid )
212     C _EXCH_XYZ_R4(recip_HFacS , myThid )
213     C _EXCH_XYZ_R4(maskW , myThid )
214     C _EXCH_XYZ_R4(maskS , myThid )
215 adcroft 1.1
216     C Calculate recipricols grid lengths
217     DO bj = myByLo(myThid), myByHi(myThid)
218     DO bi = myBxLo(myThid), myBxHi(myThid)
219 adcroft 1.19 DO J=1-Oly,sNy+Oly
220     DO I=1-Olx,sNx+Olx
221     IF ( dxG(I,J,bi,bj) .NE. 0. )
222     & recip_dxG(I,J,bi,bj)=1.d0/dxG(I,J,bi,bj)
223     IF ( dyG(I,J,bi,bj) .NE. 0. )
224     & recip_dyG(I,J,bi,bj)=1.d0/dyG(I,J,bi,bj)
225     IF ( dxC(I,J,bi,bj) .NE. 0. )
226     & recip_dxC(I,J,bi,bj)=1.d0/dxC(I,J,bi,bj)
227     IF ( dyC(I,J,bi,bj) .NE. 0. )
228     & recip_dyC(I,J,bi,bj)=1.d0/dyC(I,J,bi,bj)
229     IF ( dxF(I,J,bi,bj) .NE. 0. )
230     & recip_dxF(I,J,bi,bj)=1.d0/dxF(I,J,bi,bj)
231     IF ( dyF(I,J,bi,bj) .NE. 0. )
232     & recip_dyF(I,J,bi,bj)=1.d0/dyF(I,J,bi,bj)
233     IF ( dxV(I,J,bi,bj) .NE. 0. )
234     & recip_dxV(I,J,bi,bj)=1.d0/dxV(I,J,bi,bj)
235     IF ( dyU(I,J,bi,bj) .NE. 0. )
236     & recip_dyU(I,J,bi,bj)=1.d0/dyU(I,J,bi,bj)
237     IF ( rA(I,J,bi,bj) .NE. 0. )
238     & recip_rA(I,J,bi,bj)=1.d0/rA(I,J,bi,bj)
239     IF ( rAs(I,J,bi,bj) .NE. 0. )
240     & recip_rAs(I,J,bi,bj)=1.d0/rAs(I,J,bi,bj)
241     IF ( rAw(I,J,bi,bj) .NE. 0. )
242     & recip_rAw(I,J,bi,bj)=1.d0/rAw(I,J,bi,bj)
243     IF ( rAz(I,J,bi,bj) .NE. 0. )
244     & recip_rAz(I,J,bi,bj)=1.d0/rAz(I,J,bi,bj)
245 adcroft 1.1 ENDDO
246     ENDDO
247     ENDDO
248     ENDDO
249 adcroft 1.19 C Do not need these since above denominators are valid over full range
250     C _EXCH_XY_R4(recip_dxG, myThid )
251     C _EXCH_XY_R4(recip_dyG, myThid )
252     C _EXCH_XY_R4(recip_dxC, myThid )
253     C _EXCH_XY_R4(recip_dyC, myThid )
254     C _EXCH_XY_R4(recip_dxF, myThid )
255     C _EXCH_XY_R4(recip_dyF, myThid )
256     C _EXCH_XY_R4(recip_dxV, myThid )
257     C _EXCH_XY_R4(recip_dyU, myThid )
258     C _EXCH_XY_R4(recip_rAw, myThid )
259     C _EXCH_XY_R4(recip_rAs, myThid )
260 adcroft 1.1
261 adcroft 1.15 #ifdef ALLOW_NONHYDROSTATIC
262     C-- Calculate the reciprocal hfac distance/volume for W cells
263     DO bj = myByLo(myThid), myByHi(myThid)
264     DO bi = myBxLo(myThid), myBxHi(myThid)
265     DO K=1,Nr
266     Km1=max(K-1,1)
267     hFacUpper=drF(Km1)/(drF(Km1)+drF(K))
268     IF (Km1.EQ.K) hFacUpper=0.
269     hFacLower=drF(K)/(drF(Km1)+drF(K))
270 adcroft 1.19 DO J=1-Oly,sNy+Oly
271     DO I=1-Olx,sNx+Olx
272 adcroft 1.15 IF (hFacC(I,J,K,bi,bj).NE.0.) THEN
273     IF (hFacC(I,J,K,bi,bj).LE.0.5) THEN
274     recip_hFacU(I,J,K,bi,bj)=
275     & hFacUpper+hFacLower*hFacC(I,J,K,bi,bj)
276     ELSE
277     recip_hFacU(I,J,K,bi,bj)=1.
278     ENDIF
279     ELSE
280     recip_hFacU(I,J,K,bi,bj)=0.
281     ENDIF
282     IF (recip_hFacU(I,J,K,bi,bj).NE.0.)
283     & recip_hFacU(I,J,K,bi,bj)=1./recip_hFacU(I,J,K,bi,bj)
284     ENDDO
285     ENDDO
286     ENDDO
287     ENDDO
288     ENDDO
289 adcroft 1.19 C _EXCH_XY_R4(recip_hFacU, myThid )
290 adcroft 1.15 #endif
291 adcroft 1.1 C
292     RETURN
293     END

  ViewVC Help
Powered by ViewVC 1.1.22