/[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.31 - (show annotations) (download)
Tue Feb 7 11:47:48 2006 UTC (18 years, 3 months ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint58b_post, checkpoint58f_post, checkpoint58d_post, checkpoint58a_post, checkpoint58m_post, checkpoint58e_post, checkpoint58k_post, checkpoint58l_post, checkpoint58g_post, checkpoint58h_post, checkpoint58j_post, checkpoint58i_post, checkpoint58c_post
Changes since 1.30: +33 -1 lines
o add hooks for new package shelfice, painless

1 C $Header: /u/gcmpack/MITgcm/model/src/ini_masks_etc.F,v 1.30 2005/07/20 22:31:47 jmc Exp $
2 C $Name: $
3
4 #include "PACKAGES_CONFIG.h"
5 #include "CPP_OPTIONS.h"
6
7 CBOP
8 C !ROUTINE: INI_MASKS_ETC
9 C !INTERFACE:
10 SUBROUTINE INI_MASKS_ETC( myThid )
11 C !DESCRIPTION: \bv
12 C *==========================================================*
13 C | SUBROUTINE INI_MASKS_ETC
14 C | o Initialise masks and topography factors
15 C *==========================================================*
16 C | These arrays are used throughout the code and describe
17 C | the topography of the domain through masks (0s and 1s)
18 C | and fractional height factors (0<hFac<1). The latter
19 C | distinguish between the lopped-cell and full-step
20 C | topographic representations.
21 C *==========================================================*
22 C \ev
23
24 C !USES:
25 IMPLICIT NONE
26 C === Global variables ===
27 #include "SIZE.h"
28 #include "EEPARAMS.h"
29 #include "PARAMS.h"
30 #include "GRID.h"
31 #include "SURFACE.h"
32 #ifdef ALLOW_SHELFICE
33 # include "SHELFICE.h"
34 #endif /* ALLOW_SHELFICE */
35
36 C !INPUT/OUTPUT PARAMETERS:
37 C == Routine arguments ==
38 C myThid - Number of this instance of INI_MASKS_ETC
39 INTEGER myThid
40
41 C !LOCAL VARIABLES:
42 C == Local variables in common ==
43 C tmpfld - Temporary array used to compute & write Total Depth
44 C has to be in common for multi threading
45 COMMON / LOCAL_INI_MASKS_ETC / tmpfld
46 _RS tmpfld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
47 C == Local variables ==
48 C bi,bj - Loop counters
49 C I,J,K
50 INTEGER bi, bj
51 INTEGER I, J, K
52 #ifdef ALLOW_NONHYDROSTATIC
53 INTEGER Km1
54 _RL hFacUpper,hFacLower
55 #endif
56 _RL hFacCtmp
57 _RL hFacMnSz
58 _RL tileArea
59 CEOP
60
61 C- Calculate lopping factor hFacC : over-estimate the part inside of the domain
62 C taking into account the lower_R Boundary (Bathymetrie / Top of Atmos)
63 DO bj=myByLo(myThid), myByHi(myThid)
64 DO bi=myBxLo(myThid), myBxHi(myThid)
65 DO K=1, Nr
66 hFacMnSz=max( hFacMin, min(hFacMinDr*recip_drF(k),1. _d 0) )
67 DO J=1-Oly,sNy+Oly
68 DO I=1-Olx,sNx+Olx
69 C o Non-dimensional distance between grid bound. and domain lower_R bound.
70 hFacCtmp = (rF(K)-R_low(I,J,bi,bj))*recip_drF(K)
71 C o Select between, closed, open or partial (0,1,0-1)
72 hFacCtmp=min( max( hFacCtmp, 0. _d 0) , 1. _d 0)
73 C o Impose minimum fraction and/or size (dimensional)
74 IF (hFacCtmp.LT.hFacMnSz) THEN
75 IF (hFacCtmp.LT.hFacMnSz*0.5) THEN
76 hFacC(I,J,K,bi,bj)=0.
77 ELSE
78 hFacC(I,J,K,bi,bj)=hFacMnSz
79 ENDIF
80 ELSE
81 hFacC(I,J,K,bi,bj)=hFacCtmp
82 ENDIF
83 ENDDO
84 ENDDO
85 ENDDO
86
87 C- Re-calculate lower-R Boundary position, taking into account hFacC
88 DO J=1-Oly,sNy+Oly
89 DO I=1-Olx,sNx+Olx
90 R_low(I,J,bi,bj) = rF(1)
91 DO K=Nr,1,-1
92 R_low(I,J,bi,bj) = R_low(I,J,bi,bj)
93 & - drF(k)*hFacC(I,J,K,bi,bj)
94 ENDDO
95 ENDDO
96 ENDDO
97 C - end bi,bj loops.
98 ENDDO
99 ENDDO
100
101 C- Calculate lopping factor hFacC : Remove part outside of the domain
102 C taking into account the Reference (=at rest) Surface Position Ro_surf
103 DO bj=myByLo(myThid), myByHi(myThid)
104 DO bi=myBxLo(myThid), myBxHi(myThid)
105 DO K=1, Nr
106 hFacMnSz=max( hFacMin, min(hFacMinDr*recip_drF(k),1. _d 0) )
107 DO J=1-Oly,sNy+Oly
108 DO I=1-Olx,sNx+Olx
109 C o Non-dimensional distance between grid boundary and model surface
110 hFacCtmp = (rF(k)-Ro_surf(I,J,bi,bj))*recip_drF(K)
111 C o Reduce the previous fraction : substract the outside part.
112 hFacCtmp = hFacC(I,J,K,bi,bj) - max( hFacCtmp, 0. _d 0)
113 C o set to zero if empty Column :
114 hFacCtmp = max( hFacCtmp, 0. _d 0)
115 C o Impose minimum fraction and/or size (dimensional)
116 IF (hFacCtmp.LT.hFacMnSz) THEN
117 IF (hFacCtmp.LT.hFacMnSz*0.5) THEN
118 hFacC(I,J,K,bi,bj)=0.
119 ELSE
120 hFacC(I,J,K,bi,bj)=hFacMnSz
121 ENDIF
122 ELSE
123 hFacC(I,J,K,bi,bj)=hFacCtmp
124 ENDIF
125 ENDDO
126 ENDDO
127 ENDDO
128
129 #ifdef ALLOW_SHELFICE
130 C-- compute contributions of shelf ice to looping factors
131 IF ( useShelfIce ) THEN
132 DO K=1, Nr
133 hFacMnSz=max( hFacMin, min(hFacMinDr*recip_drF(k),1. _d 0) )
134 DO J=1-Oly,sNy+Oly
135 DO I=1-Olx,sNx+Olx
136 C o Non-dimensional distance between grid boundary and model surface
137 hFacCtmp = (rF(k)-R_shelfIce(I,J,bi,bj))*recip_drF(K)
138 C o Reduce the previous fraction : substract the outside part.
139 hFacCtmp = hFacC(I,J,K,bi,bj) - max( hFacCtmp, 0. _d 0)
140 C o set to zero if empty Column :
141 hFacCtmp = max( hFacCtmp, 0. _d 0)
142 C o Impose minimum fraction and/or size (dimensional)
143 IF (hFacCtmp.LT.hFacMnSz) THEN
144 IF (hFacCtmp.LT.hFacMnSz*0.5) THEN
145 hFacC(I,J,K,bi,bj)=0.
146 ELSE
147 hFacC(I,J,K,bi,bj)=hFacMnSz
148 ENDIF
149 ELSE
150 hFacC(I,J,K,bi,bj)=hFacCtmp
151 ENDIF
152 ENDDO
153 ENDDO
154 ENDDO
155 ENDIF
156 #endif /* ALLOW_SHELFICE */
157
158 C- Re-calculate Reference surface position, taking into account hFacC
159 C initialize Total column fluid thickness and surface k index
160 C Note: if no fluid (continent) ==> ksurf = Nr+1
161 DO J=1-Oly,sNy+Oly
162 DO I=1-Olx,sNx+Olx
163 tmpfld(I,J,bi,bj) = 0.
164 ksurfC(I,J,bi,bj) = Nr+1
165 maskH(i,j,bi,bj) = 0.
166 Ro_surf(I,J,bi,bj) = R_low(I,J,bi,bj)
167 DO K=Nr,1,-1
168 Ro_surf(I,J,bi,bj) = Ro_surf(I,J,bi,bj)
169 & + drF(k)*hFacC(I,J,K,bi,bj)
170 IF (hFacC(I,J,K,bi,bj).NE.0.) THEN
171 ksurfC(I,J,bi,bj) = k
172 maskH(i,j,bi,bj) = 1.
173 tmpfld(i,j,bi,bj) = tmpfld(i,j,bi,bj) + 1.
174 ENDIF
175 ENDDO
176 kLowC(I,J,bi,bj) = 0
177 DO K= 1, Nr
178 IF (hFacC(I,J,K,bi,bj).NE.0) THEN
179 kLowC(I,J,bi,bj) = K
180 ENDIF
181 ENDDO
182 ENDDO
183 ENDDO
184 C - end bi,bj loops.
185 ENDDO
186 ENDDO
187
188 C CALL PLOT_FIELD_XYRS( tmpfld,
189 C & 'Model Depths K Index' , 1, myThid )
190 CALL PLOT_FIELD_XYRS(R_low,
191 & 'Model R_low (ini_masks_etc)', 1, myThid)
192 CALL PLOT_FIELD_XYRS(Ro_surf,
193 & 'Model Ro_surf (ini_masks_etc)', 1, myThid)
194
195 C Calculate quantities derived from XY depth map
196 globalArea = 0. _d 0
197 DO bj = myByLo(myThid), myByHi(myThid)
198 DO bi = myBxLo(myThid), myBxHi(myThid)
199 DO j=1-Oly,sNy+Oly
200 DO i=1-Olx,sNx+Olx
201 C Total fluid column thickness (r_unit) :
202 c Rcolumn(i,j,bi,bj)= Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
203 tmpfld(i,j,bi,bj) = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
204 C Inverse of fluid column thickness (1/r_unit)
205 IF ( tmpfld(i,j,bi,bj) .LE. 0. ) THEN
206 recip_Rcol(i,j,bi,bj) = 0.
207 ELSE
208 recip_Rcol(i,j,bi,bj) = 1. / tmpfld(i,j,bi,bj)
209 ENDIF
210 ENDDO
211 ENDDO
212 C- Compute the domain Area:
213 tileArea = 0. _d 0
214 DO j=1,sNy
215 DO i=1,sNx
216 tileArea = tileArea + rA(i,j,bi,bj)*maskH(i,j,bi,bj)
217 ENDDO
218 ENDDO
219 globalArea = globalArea + tileArea
220 ENDDO
221 ENDDO
222 C _EXCH_XY_R4( recip_Rcol, myThid )
223 _GLOBAL_SUM_R8( globalArea, myThid )
224
225 C hFacW and hFacS (at U and V points)
226 DO bj=myByLo(myThid), myByHi(myThid)
227 DO bi=myBxLo(myThid), myBxHi(myThid)
228 DO K=1, Nr
229 DO J=1-Oly,sNy+Oly
230 DO I=2-Olx,sNx+Olx
231 hFacW(I,J,K,bi,bj)=
232 & MIN(hFacC(I,J,K,bi,bj),hFacC(I-1,J,K,bi,bj))
233 ENDDO
234 ENDDO
235 DO J=2-Oly,sNy+oly
236 DO I=1-Olx,sNx+Olx
237 hFacS(I,J,K,bi,bj)=
238 & MIN(hFacC(I,J,K,bi,bj),hFacC(I,J-1,K,bi,bj))
239 ENDDO
240 ENDDO
241 ENDDO
242 ENDDO
243 ENDDO
244 CALL EXCH_UV_XYZ_RS(hFacW,hFacS,.FALSE.,myThid)
245 C The following block allows thin walls representation of non-periodic
246 C boundaries such as happen on the lat-lon grid at the N/S poles.
247 C We should really supply a flag for doing this.
248 DO bj=myByLo(myThid), myByHi(myThid)
249 DO bi=myBxLo(myThid), myBxHi(myThid)
250 DO K=1, Nr
251 DO J=1-Oly,sNy+Oly
252 DO I=1-Olx,sNx+Olx
253 IF (DYG(I,J,bi,bj).EQ.0.) hFacW(I,J,K,bi,bj)=0.
254 IF (DXG(I,J,bi,bj).EQ.0.) hFacS(I,J,K,bi,bj)=0.
255 ENDDO
256 ENDDO
257 ENDDO
258 ENDDO
259 ENDDO
260
261 C- Write to disk: Total Column Thickness & hFac(C,W,S):
262 _BARRIER
263 c _BEGIN_MASTER( myThid )
264 C This I/O is now done in write_grid.F
265 C CALL MDSWRITEFIELD( 'Depth', writeBinaryPrec, .TRUE.,
266 C & 'RS', 1, tmpfld, 1, -1, myThid )
267 c CALL WRITE_FLD_XY_RS( 'Depth',' ',tmpfld,0,myThid)
268 c CALL WRITE_FLD_XYZ_RS( 'hFacC',' ',hFacC,0,myThid)
269 c CALL WRITE_FLD_XYZ_RS( 'hFacW',' ',hFacW,0,myThid)
270 c CALL WRITE_FLD_XYZ_RS( 'hFacS',' ',hFacS,0,myThid)
271 c _END_MASTER(myThid)
272
273 CALL PLOT_FIELD_XYZRS( hFacC, 'hFacC' , Nr, 1, myThid )
274 CALL PLOT_FIELD_XYZRS( hFacW, 'hFacW' , Nr, 1, myThid )
275 CALL PLOT_FIELD_XYZRS( hFacS, 'hFacS' , Nr, 1, myThid )
276
277 C Masks and reciprocals of hFac[CWS]
278 DO bj = myByLo(myThid), myByHi(myThid)
279 DO bi = myBxLo(myThid), myBxHi(myThid)
280 DO K=1,Nr
281 DO J=1-Oly,sNy+Oly
282 DO I=1-Olx,sNx+Olx
283 IF (HFacC(I,J,K,bi,bj) .NE. 0. ) THEN
284 recip_HFacC(I,J,K,bi,bj) = 1. / HFacC(I,J,K,bi,bj)
285 maskC(I,J,K,bi,bj) = 1.
286 ELSE
287 recip_HFacC(I,J,K,bi,bj) = 0.
288 maskC(I,J,K,bi,bj) = 0.
289 ENDIF
290 IF (HFacW(I,J,K,bi,bj) .NE. 0. ) THEN
291 recip_HFacW(I,J,K,bi,bj) = 1. / HFacW(I,J,K,bi,bj)
292 maskW(I,J,K,bi,bj) = 1.
293 ELSE
294 recip_HFacW(I,J,K,bi,bj) = 0.
295 maskW(I,J,K,bi,bj) = 0.
296 ENDIF
297 IF (HFacS(I,J,K,bi,bj) .NE. 0. ) THEN
298 recip_HFacS(I,J,K,bi,bj) = 1. / HFacS(I,J,K,bi,bj)
299 maskS(I,J,K,bi,bj) = 1.
300 ELSE
301 recip_HFacS(I,J,K,bi,bj) = 0.
302 maskS(I,J,K,bi,bj) = 0.
303 ENDIF
304 ENDDO
305 ENDDO
306 ENDDO
307 C- Calculate surface k index for interface W & S (U & V points)
308 DO J=1-Oly,sNy+Oly
309 DO I=1-Olx,sNx+Olx
310 ksurfW(I,J,bi,bj) = Nr+1
311 ksurfS(I,J,bi,bj) = Nr+1
312 DO k=Nr,1,-1
313 IF (hFacW(I,J,K,bi,bj).NE.0.) ksurfW(I,J,bi,bj) = k
314 IF (hFacS(I,J,K,bi,bj).NE.0.) ksurfS(I,J,bi,bj) = k
315 ENDDO
316 ENDDO
317 ENDDO
318 C - end bi,bj loops.
319 ENDDO
320 ENDDO
321 C _EXCH_XYZ_R4(recip_HFacC , myThid )
322 C _EXCH_XYZ_R4(recip_HFacW , myThid )
323 C _EXCH_XYZ_R4(recip_HFacS , myThid )
324 C _EXCH_XYZ_R4(maskW , myThid )
325 C _EXCH_XYZ_R4(maskS , myThid )
326
327 C Calculate recipricols grid lengths
328 DO bj = myByLo(myThid), myByHi(myThid)
329 DO bi = myBxLo(myThid), myBxHi(myThid)
330 DO J=1-Oly,sNy+Oly
331 DO I=1-Olx,sNx+Olx
332 IF ( dxG(I,J,bi,bj) .NE. 0. )
333 & recip_dxG(I,J,bi,bj)=1.d0/dxG(I,J,bi,bj)
334 IF ( dyG(I,J,bi,bj) .NE. 0. )
335 & recip_dyG(I,J,bi,bj)=1.d0/dyG(I,J,bi,bj)
336 IF ( dxC(I,J,bi,bj) .NE. 0. )
337 & recip_dxC(I,J,bi,bj)=1.d0/dxC(I,J,bi,bj)
338 IF ( dyC(I,J,bi,bj) .NE. 0. )
339 & recip_dyC(I,J,bi,bj)=1.d0/dyC(I,J,bi,bj)
340 IF ( dxF(I,J,bi,bj) .NE. 0. )
341 & recip_dxF(I,J,bi,bj)=1.d0/dxF(I,J,bi,bj)
342 IF ( dyF(I,J,bi,bj) .NE. 0. )
343 & recip_dyF(I,J,bi,bj)=1.d0/dyF(I,J,bi,bj)
344 IF ( dxV(I,J,bi,bj) .NE. 0. )
345 & recip_dxV(I,J,bi,bj)=1.d0/dxV(I,J,bi,bj)
346 IF ( dyU(I,J,bi,bj) .NE. 0. )
347 & recip_dyU(I,J,bi,bj)=1.d0/dyU(I,J,bi,bj)
348 IF ( rA(I,J,bi,bj) .NE. 0. )
349 & recip_rA(I,J,bi,bj)=1.d0/rA(I,J,bi,bj)
350 IF ( rAs(I,J,bi,bj) .NE. 0. )
351 & recip_rAs(I,J,bi,bj)=1.d0/rAs(I,J,bi,bj)
352 IF ( rAw(I,J,bi,bj) .NE. 0. )
353 & recip_rAw(I,J,bi,bj)=1.d0/rAw(I,J,bi,bj)
354 IF ( rAz(I,J,bi,bj) .NE. 0. )
355 & recip_rAz(I,J,bi,bj)=1.d0/rAz(I,J,bi,bj)
356 ENDDO
357 ENDDO
358 ENDDO
359 ENDDO
360 C Do not need these since above denominators are valid over full range
361 C _EXCH_XY_R4(recip_dxG, myThid )
362 C _EXCH_XY_R4(recip_dyG, myThid )
363 C _EXCH_XY_R4(recip_dxC, myThid )
364 C _EXCH_XY_R4(recip_dyC, myThid )
365 C _EXCH_XY_R4(recip_dxF, myThid )
366 C _EXCH_XY_R4(recip_dyF, myThid )
367 C _EXCH_XY_R4(recip_dxV, myThid )
368 C _EXCH_XY_R4(recip_dyU, myThid )
369 C _EXCH_XY_R4(recip_rAw, myThid )
370 C _EXCH_XY_R4(recip_rAs, myThid )
371
372 #ifdef ALLOW_NONHYDROSTATIC
373 C-- Calculate the reciprocal hfac distance/volume for W cells
374 DO bj = myByLo(myThid), myByHi(myThid)
375 DO bi = myBxLo(myThid), myBxHi(myThid)
376 DO K=1,Nr
377 Km1=max(K-1,1)
378 hFacUpper=drF(Km1)/(drF(Km1)+drF(K))
379 IF (Km1.EQ.K) hFacUpper=0.
380 hFacLower=drF(K)/(drF(Km1)+drF(K))
381 DO J=1-Oly,sNy+Oly
382 DO I=1-Olx,sNx+Olx
383 IF (hFacC(I,J,K,bi,bj).NE.0.) THEN
384 IF (hFacC(I,J,K,bi,bj).LE.0.5) THEN
385 recip_hFacU(I,J,K,bi,bj)=
386 & hFacUpper+hFacLower*hFacC(I,J,K,bi,bj)
387 ELSE
388 recip_hFacU(I,J,K,bi,bj)=1.
389 ENDIF
390 ELSE
391 recip_hFacU(I,J,K,bi,bj)=0.
392 ENDIF
393 IF (recip_hFacU(I,J,K,bi,bj).NE.0.)
394 & recip_hFacU(I,J,K,bi,bj)=1./recip_hFacU(I,J,K,bi,bj)
395 ENDDO
396 ENDDO
397 ENDDO
398 ENDDO
399 ENDDO
400 C _EXCH_XY_R4(recip_hFacU, myThid )
401 #endif
402 C
403 RETURN
404 END

  ViewVC Help
Powered by ViewVC 1.1.22