/[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.23 - (hide annotations) (download)
Mon Aug 27 18:42:37 2001 UTC (22 years, 8 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint40pre9, checkpoint40
Changes since 1.22: +16 -3 lines
Change name and definition of k_surf (now ksurfC)
  compute also ksurfW,ksurfS (for NonLin-FreeSurf)

1 jmc 1.23 C $Header: /u/gcmpack/models/MITgcmUV/model/src/ini_masks_etc.F,v 1.22 2001/07/05 21:44:25 jmc Exp $
2 adcroft 1.21 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 adcroft 1.21 #include "SURFACE.h"
26 adcroft 1.1
27     C == Routine arguments ==
28 cnh 1.6 C myThid - Number of this instance of INI_MASKS_ETC
29 adcroft 1.1 INTEGER myThid
30     CEndOfInterface
31    
32 jmc 1.22 C == Local variables in common ==
33     C tmpfld - Temporary array used to compute & write Total Depth
34     C has to be in common for multi threading
35     COMMON / LOCAL_INI_MASKS_ETC / tmpfld
36     _RS tmpfld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
37    
38 adcroft 1.1 C == Local variables ==
39     C bi,bj - Loop counters
40     C I,J,K
41     INTEGER bi, bj
42     INTEGER I, J, K
43 adcroft 1.15 #ifdef ALLOW_NONHYDROSTATIC
44     INTEGER Km1
45     _RL hFacUpper,hFacLower
46     #endif
47 adcroft 1.21 _RL hFacCtmp
48 adcroft 1.19 _RL hFacMnSz
49    
50 adcroft 1.21 C- Calculate lopping factor hFacC : over-estimate the part inside of the domain
51     C taking into account the lower_R Boundary (Bathymetrie / Top of Atmos)
52 adcroft 1.2 DO bj=myByLo(myThid), myByHi(myThid)
53     DO bi=myBxLo(myThid), myBxHi(myThid)
54 cnh 1.4 DO K=1, Nr
55 adcroft 1.21 hFacMnSz=max( hFacMin, min(hFacMinDr*recip_drF(k),1. _d 0) )
56 adcroft 1.19 DO J=1-Oly,sNy+Oly
57     DO I=1-Olx,sNx+Olx
58 adcroft 1.21 C o Non-dimensional distance between grid bound. and domain lower_R bound.
59     hFacCtmp = (rF(K)-R_low(I,J,bi,bj))*recip_drF(K)
60     C o Select between, closed, open or partial (0,1,0-1)
61 adcroft 1.19 hFacCtmp=min( max( hFacCtmp, 0. _d 0) , 1. _d 0)
62 adcroft 1.21 C o Impose minimum fraction and/or size (dimensional)
63     IF (hFacCtmp.LT.hFacMnSz) THEN
64     IF (hFacCtmp.LT.hFacMnSz*0.5) THEN
65 adcroft 1.3 hFacC(I,J,K,bi,bj)=0.
66     ELSE
67 adcroft 1.19 hFacC(I,J,K,bi,bj)=hFacMnSz
68 adcroft 1.3 ENDIF
69 adcroft 1.21 ELSE
70     hFacC(I,J,K,bi,bj)=hFacCtmp
71 adcroft 1.2 ENDIF
72     ENDDO
73     ENDDO
74     ENDDO
75 adcroft 1.21
76     C- Re-calculate lower-R Boundary position, taking into account hFacC
77     DO J=1-Oly,sNy+Oly
78     DO I=1-Olx,sNx+Olx
79     R_low(I,J,bi,bj) = rF(1)
80     DO K=Nr,1,-1
81     R_low(I,J,bi,bj) = R_low(I,J,bi,bj)
82     & - drF(k)*hFacC(I,J,K,bi,bj)
83     ENDDO
84     ENDDO
85     ENDDO
86     C - end bi,bj loops.
87 adcroft 1.2 ENDDO
88     ENDDO
89 cnh 1.7
90 adcroft 1.21 C- Calculate lopping factor hFacC : Remove part outside of the domain
91     C taking into account the Reference (=at rest) Surface Position Ro_surf
92 adcroft 1.16 DO bj=myByLo(myThid), myByHi(myThid)
93     DO bi=myBxLo(myThid), myBxHi(myThid)
94 adcroft 1.21 DO K=1, Nr
95     hFacMnSz=max( hFacMin, min(hFacMinDr*recip_drF(k),1. _d 0) )
96     DO J=1-Oly,sNy+Oly
97     DO I=1-Olx,sNx+Olx
98     C o Non-dimensional distance between grid boundary and model surface
99     hFacCtmp = (rF(k)-Ro_surf(I,J,bi,bj))*recip_drF(K)
100     C o Reduce the previous fraction : substract the outside part.
101     hFacCtmp = hFacC(I,J,K,bi,bj) - max( hFacCtmp, 0. _d 0)
102     C o set to zero if empty Column :
103     hFacCtmp = max( hFacCtmp, 0. _d 0)
104     C o Impose minimum fraction and/or size (dimensional)
105     IF (hFacCtmp.LT.hFacMnSz) THEN
106     IF (hFacCtmp.LT.hFacMnSz*0.5) THEN
107     hFacC(I,J,K,bi,bj)=0.
108     ELSE
109     hFacC(I,J,K,bi,bj)=hFacMnSz
110     ENDIF
111     ELSE
112     hFacC(I,J,K,bi,bj)=hFacCtmp
113     ENDIF
114     ENDDO
115     ENDDO
116     ENDDO
117    
118     C- Re-calculate Reference surface position, taking into account hFacC
119     C initialize Total column fluid thickness and surface k index
120 jmc 1.23 C Note: if no fluid (continent) ==> ksurf = Nr+1
121 adcroft 1.19 DO J=1-Oly,sNy+Oly
122     DO I=1-Olx,sNx+Olx
123 adcroft 1.21 tmpfld(I,J,bi,bj) = 0.
124 jmc 1.23 ksurfC(I,J,bi,bj) = Nr+1
125 adcroft 1.21 Ro_surf(I,J,bi,bj) = R_low(I,J,bi,bj)
126     DO K=Nr,1,-1
127     Ro_surf(I,J,bi,bj) = Ro_surf(I,J,bi,bj)
128     & + drF(k)*hFacC(I,J,K,bi,bj)
129     IF (hFacC(I,J,K,bi,bj).NE.0.) THEN
130 jmc 1.23 ksurfC(I,J,bi,bj) = k
131 adcroft 1.21 tmpfld(i,j,bi,bj) = tmpfld(i,j,bi,bj) + 1.
132     ENDIF
133 adcroft 1.16 ENDDO
134     ENDDO
135     ENDDO
136 adcroft 1.21 C - end bi,bj loops.
137 adcroft 1.16 ENDDO
138     ENDDO
139 adcroft 1.21
140     C CALL PLOT_FIELD_XYRS( tmpfld,
141     C & 'Model Depths K Index' , 1, myThid )
142     CALL PLOT_FIELD_XYRS(R_low,
143     & 'Model R_low (ini_masks_etc)', 1, myThid)
144     CALL PLOT_FIELD_XYRS(Ro_surf,
145     & 'Model Ro_surf (ini_masks_etc)', 1, myThid)
146 adcroft 1.16
147     C Calculate quantities derived from XY depth map
148     DO bj = myByLo(myThid), myByHi(myThid)
149     DO bi = myBxLo(myThid), myBxHi(myThid)
150 adcroft 1.21 DO j=1-Oly,sNy+Oly
151     DO i=1-Olx,sNx+Olx
152     C Total fluid column thickness (r_unit) :
153     c Rcolumn(i,j,bi,bj)= Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
154     tmpfld(i,j,bi,bj) = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
155     C Inverse of fluid column thickness (1/r_unit)
156     IF ( tmpfld(i,j,bi,bj) .LE. 0. ) THEN
157     recip_Rcol(i,j,bi,bj) = 0.
158 adcroft 1.16 ELSE
159 adcroft 1.21 recip_Rcol(i,j,bi,bj) = 1. / tmpfld(i,j,bi,bj)
160 adcroft 1.16 ENDIF
161     ENDDO
162     ENDDO
163     ENDDO
164     ENDDO
165 adcroft 1.21 C _EXCH_XY_R4( recip_Rcol, myThid )
166 adcroft 1.1
167     C hFacW and hFacS (at U and V points)
168     DO bj=myByLo(myThid), myByHi(myThid)
169     DO bi=myBxLo(myThid), myBxHi(myThid)
170 cnh 1.4 DO K=1, Nr
171 adcroft 1.1 DO J=1,sNy
172     DO I=1,sNx
173     hFacW(I,J,K,bi,bj)=
174     & MIN(hFacC(I,J,K,bi,bj),hFacC(I-1,J,K,bi,bj))
175     hFacS(I,J,K,bi,bj)=
176     & MIN(hFacC(I,J,K,bi,bj),hFacC(I,J-1,K,bi,bj))
177     ENDDO
178     ENDDO
179     ENDDO
180     ENDDO
181     ENDDO
182 adcroft 1.21 CALL EXCH_UV_XYZ_RS(hFacW,hFacS,.FALSE.,myThid)
183     C The following block allows thin walls representation of non-periodic
184     C boundaries such as happen on the lat-lon grid at the N/S poles.
185     C We should really supply a flag for doing this.
186 adcroft 1.19 DO bj=myByLo(myThid), myByHi(myThid)
187     DO bi=myBxLo(myThid), myBxHi(myThid)
188     DO K=1, Nr
189     DO J=1-Oly,sNy+Oly
190 adcroft 1.21 DO I=1-Olx,sNx+Olx
191     IF (DYG(I,J,bi,bj).EQ.0.) hFacW(I,J,K,bi,bj)=0.
192     IF (DXG(I,J,bi,bj).EQ.0.) hFacS(I,J,K,bi,bj)=0.
193 adcroft 1.19 ENDDO
194     ENDDO
195     ENDDO
196     ENDDO
197     ENDDO
198 jmc 1.22
199     C- Write to disk: Total Column Thickness & hFac(C,W,S):
200     _BARRIER
201     _BEGIN_MASTER( myThid )
202     C CALL MDSWRITEFIELD( 'Depth', writeBinaryPrec, .TRUE.,
203     C & 'RS', 1, tmpfld, 1, -1, myThid )
204     CALL WRITE_FLD_XY_RS( 'Depth',' ',tmpfld,0,myThid)
205     CALL WRITE_FLD_XYZ_RS( 'hFacC',' ',hFacC,0,myThid)
206     CALL WRITE_FLD_XYZ_RS( 'hFacW',' ',hFacW,0,myThid)
207     CALL WRITE_FLD_XYZ_RS( 'hFacS',' ',hFacS,0,myThid)
208     _END_MASTER(myThid)
209 adcroft 1.19
210     CALL PLOT_FIELD_XYZRS( hFacC, 'hFacC' , Nr, 1, myThid )
211     CALL PLOT_FIELD_XYZRS( hFacW, 'hFacW' , Nr, 1, myThid )
212     CALL PLOT_FIELD_XYZRS( hFacS, 'hFacS' , Nr, 1, myThid )
213 adcroft 1.1
214     C Masks and reciprocals of hFac[CWS]
215     DO bj = myByLo(myThid), myByHi(myThid)
216     DO bi = myBxLo(myThid), myBxHi(myThid)
217 cnh 1.4 DO K=1,Nr
218 adcroft 1.19 DO J=1-Oly,sNy+Oly
219     DO I=1-Olx,sNx+Olx
220     IF (HFacC(I,J,K,bi,bj) .NE. 0. ) THEN
221     recip_HFacC(I,J,K,bi,bj) = 1. / HFacC(I,J,K,bi,bj)
222 adcroft 1.21 maskC(I,J,K,bi,bj) = 1.
223 adcroft 1.1 ELSE
224 adcroft 1.19 recip_HFacC(I,J,K,bi,bj) = 0.
225 adcroft 1.21 maskC(I,J,K,bi,bj) = 0.
226 adcroft 1.1 ENDIF
227 adcroft 1.19 IF (HFacW(I,J,K,bi,bj) .NE. 0. ) THEN
228     recip_HFacW(I,J,K,bi,bj) = 1. / HFacW(I,J,K,bi,bj)
229     maskW(I,J,K,bi,bj) = 1.
230 adcroft 1.1 ELSE
231 adcroft 1.19 recip_HFacW(I,J,K,bi,bj) = 0.
232     maskW(I,J,K,bi,bj) = 0.
233 adcroft 1.1 ENDIF
234 adcroft 1.19 IF (HFacS(I,J,K,bi,bj) .NE. 0. ) THEN
235     recip_HFacS(I,J,K,bi,bj) = 1. / HFacS(I,J,K,bi,bj)
236     maskS(I,J,K,bi,bj) = 1.
237 adcroft 1.1 ELSE
238 adcroft 1.19 recip_HFacS(I,J,K,bi,bj) = 0.
239     maskS(I,J,K,bi,bj) = 0.
240 adcroft 1.1 ENDIF
241     ENDDO
242     ENDDO
243     ENDDO
244 jmc 1.23 C- Calculate surface k index for interface W & S (U & V points)
245     DO J=1-Oly,sNy+Oly
246     DO I=1-Olx,sNx+Olx
247     ksurfW(I,J,bi,bj) = Nr+1
248     ksurfS(I,J,bi,bj) = Nr+1
249     DO k=Nr,1,-1
250     IF (hFacW(I,J,K,bi,bj).NE.0.) ksurfW(I,J,bi,bj) = k
251     IF (hFacS(I,J,K,bi,bj).NE.0.) ksurfS(I,J,bi,bj) = k
252     ENDDO
253     ENDDO
254     ENDDO
255     C - end bi,bj loops.
256 adcroft 1.1 ENDDO
257     ENDDO
258 adcroft 1.19 C _EXCH_XYZ_R4(recip_HFacC , myThid )
259     C _EXCH_XYZ_R4(recip_HFacW , myThid )
260     C _EXCH_XYZ_R4(recip_HFacS , myThid )
261     C _EXCH_XYZ_R4(maskW , myThid )
262     C _EXCH_XYZ_R4(maskS , myThid )
263 adcroft 1.1
264     C Calculate recipricols grid lengths
265     DO bj = myByLo(myThid), myByHi(myThid)
266     DO bi = myBxLo(myThid), myBxHi(myThid)
267 adcroft 1.19 DO J=1-Oly,sNy+Oly
268     DO I=1-Olx,sNx+Olx
269     IF ( dxG(I,J,bi,bj) .NE. 0. )
270     & recip_dxG(I,J,bi,bj)=1.d0/dxG(I,J,bi,bj)
271     IF ( dyG(I,J,bi,bj) .NE. 0. )
272     & recip_dyG(I,J,bi,bj)=1.d0/dyG(I,J,bi,bj)
273     IF ( dxC(I,J,bi,bj) .NE. 0. )
274     & recip_dxC(I,J,bi,bj)=1.d0/dxC(I,J,bi,bj)
275     IF ( dyC(I,J,bi,bj) .NE. 0. )
276     & recip_dyC(I,J,bi,bj)=1.d0/dyC(I,J,bi,bj)
277     IF ( dxF(I,J,bi,bj) .NE. 0. )
278     & recip_dxF(I,J,bi,bj)=1.d0/dxF(I,J,bi,bj)
279     IF ( dyF(I,J,bi,bj) .NE. 0. )
280     & recip_dyF(I,J,bi,bj)=1.d0/dyF(I,J,bi,bj)
281     IF ( dxV(I,J,bi,bj) .NE. 0. )
282     & recip_dxV(I,J,bi,bj)=1.d0/dxV(I,J,bi,bj)
283     IF ( dyU(I,J,bi,bj) .NE. 0. )
284     & recip_dyU(I,J,bi,bj)=1.d0/dyU(I,J,bi,bj)
285     IF ( rA(I,J,bi,bj) .NE. 0. )
286     & recip_rA(I,J,bi,bj)=1.d0/rA(I,J,bi,bj)
287     IF ( rAs(I,J,bi,bj) .NE. 0. )
288     & recip_rAs(I,J,bi,bj)=1.d0/rAs(I,J,bi,bj)
289     IF ( rAw(I,J,bi,bj) .NE. 0. )
290     & recip_rAw(I,J,bi,bj)=1.d0/rAw(I,J,bi,bj)
291     IF ( rAz(I,J,bi,bj) .NE. 0. )
292     & recip_rAz(I,J,bi,bj)=1.d0/rAz(I,J,bi,bj)
293 adcroft 1.1 ENDDO
294     ENDDO
295     ENDDO
296     ENDDO
297 adcroft 1.19 C Do not need these since above denominators are valid over full range
298     C _EXCH_XY_R4(recip_dxG, myThid )
299     C _EXCH_XY_R4(recip_dyG, myThid )
300     C _EXCH_XY_R4(recip_dxC, myThid )
301     C _EXCH_XY_R4(recip_dyC, myThid )
302     C _EXCH_XY_R4(recip_dxF, myThid )
303     C _EXCH_XY_R4(recip_dyF, myThid )
304     C _EXCH_XY_R4(recip_dxV, myThid )
305     C _EXCH_XY_R4(recip_dyU, myThid )
306     C _EXCH_XY_R4(recip_rAw, myThid )
307     C _EXCH_XY_R4(recip_rAs, myThid )
308 adcroft 1.1
309 adcroft 1.15 #ifdef ALLOW_NONHYDROSTATIC
310     C-- Calculate the reciprocal hfac distance/volume for W cells
311     DO bj = myByLo(myThid), myByHi(myThid)
312     DO bi = myBxLo(myThid), myBxHi(myThid)
313     DO K=1,Nr
314     Km1=max(K-1,1)
315     hFacUpper=drF(Km1)/(drF(Km1)+drF(K))
316     IF (Km1.EQ.K) hFacUpper=0.
317     hFacLower=drF(K)/(drF(Km1)+drF(K))
318 adcroft 1.19 DO J=1-Oly,sNy+Oly
319     DO I=1-Olx,sNx+Olx
320 adcroft 1.15 IF (hFacC(I,J,K,bi,bj).NE.0.) THEN
321     IF (hFacC(I,J,K,bi,bj).LE.0.5) THEN
322     recip_hFacU(I,J,K,bi,bj)=
323     & hFacUpper+hFacLower*hFacC(I,J,K,bi,bj)
324     ELSE
325     recip_hFacU(I,J,K,bi,bj)=1.
326     ENDIF
327     ELSE
328     recip_hFacU(I,J,K,bi,bj)=0.
329     ENDIF
330     IF (recip_hFacU(I,J,K,bi,bj).NE.0.)
331     & recip_hFacU(I,J,K,bi,bj)=1./recip_hFacU(I,J,K,bi,bj)
332     ENDDO
333     ENDDO
334     ENDDO
335     ENDDO
336     ENDDO
337 adcroft 1.19 C _EXCH_XY_R4(recip_hFacU, myThid )
338 adcroft 1.15 #endif
339 adcroft 1.1 C
340     RETURN
341     END

  ViewVC Help
Powered by ViewVC 1.1.22