/[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.19 - (hide annotations) (download)
Fri Feb 2 21:04:48 2001 UTC (23 years, 4 months ago) by adcroft
Branch: MAIN
Changes since 1.18: +139 -85 lines
Merged changes from branch "branch-atmos-merge" into MAIN (checkpoint34)
 - substantial modifications to algorithm sequence (dynamics.F)
 - packaged OBCS, Shapiro filter, Zonal filter, Atmospheric Physics

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

  ViewVC Help
Powered by ViewVC 1.1.22