/[MITgcm]/MITgcm_contrib/dgoldberg/shelfice_remeshing_old/code/ini_masks_etc.F
ViewVC logotype

Contents of /MITgcm_contrib/dgoldberg/shelfice_remeshing_old/code/ini_masks_etc.F

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


Revision 1.1 - (show annotations) (download)
Thu Jul 19 11:14:58 2018 UTC (7 years ago) by dgoldberg
Branch: MAIN
CVS Tags: HEAD
record of verification_other experiment before modification for github

1 C $Header: /u/gcmpack/MITgcm_contrib/verification_other/shelfice_remeshing/code/ini_masks_etc.F,v 1.3 2017/04/10 23:55:33 jmc Exp $
2 C $Name: $
3
4 #include "PACKAGES_CONFIG.h"
5 #include "CPP_OPTIONS.h"
6 #ifdef ALLOW_SHELFICE
7 # include "SHELFICE_OPTIONS.h"
8 #endif /* ALLOW_SHELFICE */
9
10 CBOP
11 C !ROUTINE: INI_MASKS_ETC
12 C !INTERFACE:
13 SUBROUTINE INI_MASKS_ETC( myThid )
14 C !DESCRIPTION: \bv
15 C *==========================================================*
16 C | SUBROUTINE INI_MASKS_ETC
17 C | o Initialise masks and topography factors
18 C *==========================================================*
19 C | These arrays are used throughout the code and describe
20 C | the topography of the domain through masks (0s and 1s)
21 C | and fractional height factors (0<hFac<1). The latter
22 C | distinguish between the lopped-cell and full-step
23 C | topographic representations.
24 C *==========================================================*
25 C \ev
26
27 C !USES:
28 IMPLICIT NONE
29 C === Global variables ===
30 #include "SIZE.h"
31 #include "EEPARAMS.h"
32 #include "PARAMS.h"
33 #include "GRID.h"
34 #ifdef NONLIN_FRSURF
35 # include "SURFACE.h"
36 #endif /* NONLIN_FRSURF */
37
38 C !INPUT/OUTPUT PARAMETERS:
39 C == Routine arguments ==
40 C myThid :: Number of this instance of INI_MASKS_ETC
41 INTEGER myThid
42
43 C !LOCAL VARIABLES:
44 C == Local variables ==
45 C bi,bj :: tile indices
46 C i,j,k :: Loop counters
47 C tmpFld :: Temporary array used to compute & write Total Depth
48 _RS tmpFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
49 INTEGER bi, bj
50 INTEGER i, j, k
51 _RL hFacCtmp
52 _RL hFacMnSz
53 CEOP
54
55 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
56
57 #ifdef ALLOW_SHELFICE
58 IF ( useShelfIce ) THEN
59 C-- Modify ocean upper boundary position according to ice-shelf topography
60 CALL SHELFICE_INIT_DEPTHS(
61 U R_low, Ro_surf,
62 I myThid )
63 ENDIF
64 #endif /* ALLOW_SHELFICE */
65
66 IF ( selectSigmaCoord.EQ.0 ) THEN
67 C--- r-coordinate with partial-cell or full cell mask
68
69 C-- Calculate lopping factor hFacC : over-estimate the part inside of the domain
70 C taking into account the lower_R Boundary (Bathymetrie / Top of Atmos)
71 DO bj=myByLo(myThid), myByHi(myThid)
72 DO bi=myBxLo(myThid), myBxHi(myThid)
73 DO k=1, Nr
74 hFacMnSz = MAX( hFacMin, MIN(hFacMinDr*recip_drF(k),oneRL) )
75 DO j=1-OLy,sNy+OLy
76 DO i=1-OLx,sNx+OLx
77 C o Non-dimensional distance between grid bound. and domain lower_R bound.
78 hFacCtmp = (rF(k)-R_low(i,j,bi,bj))*recip_drF(k)
79 C o Select between, closed, open or partial (0,1,0-1)
80 hFacCtmp = MIN( MAX( hFacCtmp, zeroRL ) , oneRL )
81 C o Impose minimum fraction and/or size (dimensional)
82 IF ( hFacCtmp.LT.hFacMnSz ) THEN
83 IF ( hFacCtmp.LT.hFacMnSz*0.5 ) THEN
84 hFacC(i,j,k,bi,bj) = 0.
85 ELSE
86 hFacC(i,j,k,bi,bj) = hFacMnSz
87 ENDIF
88 ELSE
89 hFacC(i,j,k,bi,bj) = hFacCtmp
90 ENDIF
91 ENDDO
92 ENDDO
93 ENDDO
94
95 C- Re-calculate lower-R Boundary position, taking into account hFacC
96 DO j=1-OLy,sNy+OLy
97 DO i=1-OLx,sNx+OLx
98 R_low(i,j,bi,bj) = rF(1)
99 ENDDO
100 ENDDO
101 DO k=Nr,1,-1
102 DO j=1-OLy,sNy+OLy
103 DO i=1-OLx,sNx+OLx
104 R_low(i,j,bi,bj) = R_low(i,j,bi,bj)
105 & - drF(k)*hFacC(i,j,k,bi,bj)
106 ENDDO
107 ENDDO
108 ENDDO
109 C- end bi,bj loops.
110 ENDDO
111 ENDDO
112
113 C-- Calculate lopping factor hFacC : Remove part outside of the domain
114 C taking into account the Reference (=at rest) Surface Position Ro_surf
115 DO bj=myByLo(myThid), myByHi(myThid)
116 DO bi=myBxLo(myThid), myBxHi(myThid)
117 DO k=1, Nr
118 hFacMnSz = MAX( hFacMin, MIN(hFacMinDr*recip_drF(k),oneRL) )
119 DO j=1-OLy,sNy+OLy
120 DO i=1-OLx,sNx+OLx
121 C o Non-dimensional distance between grid boundary and model surface
122 hFacCtmp = (rF(k)-Ro_surf(i,j,bi,bj))*recip_drF(k)
123 C o Reduce the previous fraction : substract the outside part.
124 hFacCtmp = hFacC(i,j,k,bi,bj) - MAX( hFacCtmp, zeroRL )
125 C o set to zero if empty Column :
126 hFacCtmp = MAX( hFacCtmp, zeroRL )
127 C o Impose minimum fraction and/or size (dimensional)
128 IF ( hFacCtmp.LT.hFacMnSz ) THEN
129 IF ( hFacCtmp.LT.hFacMnSz*0.5 ) THEN
130 hFacC(i,j,k,bi,bj) = 0.
131 ELSE
132 hFacC(i,j,k,bi,bj) = hFacMnSz
133 ENDIF
134 ELSE
135 hFacC(i,j,k,bi,bj) = hFacCtmp
136 ENDIF
137 ENDDO
138 ENDDO
139 ENDDO
140 ENDDO
141 ENDDO
142
143 c#ifdef ALLOW_SHELFICE
144 c IF ( useShelfIce ) THEN
145 C-- Modify lopping factor hFacC : Remove part outside of the domain
146 C taking into account the Reference (=at rest) Surface Position Ro_shelfIce
147 c CALL SHELFICE_UPDATE_MASKS(
148 c I rF, recip_drF,
149 c U hFacC,
150 c I myThid )
151 c ENDIF
152 c#endif /* ALLOW_SHELFICE */
153
154 C- Re-calculate Reference surface position, taking into account hFacC
155 C initialize Total column fluid thickness and surface k index
156 C Note: if no fluid (continent) ==> kSurf = Nr+1
157 DO bj=myByLo(myThid), myByHi(myThid)
158 DO bi=myBxLo(myThid), myBxHi(myThid)
159 DO j=1-OLy,sNy+OLy
160 DO i=1-OLx,sNx+OLx
161 tmpFld(i,j,bi,bj) = 0.
162 kSurfC(i,j,bi,bj) = Nr+1
163 Ro_surf(i,j,bi,bj) = R_low(i,j,bi,bj)
164 DO k=Nr,1,-1
165 Ro_surf(i,j,bi,bj) = Ro_surf(i,j,bi,bj)
166 & + drF(k)*hFacC(i,j,k,bi,bj)
167 IF ( hFacC(i,j,k,bi,bj).NE.zeroRS ) THEN
168 kSurfC(i,j,bi,bj) = k
169 tmpFld(i,j,bi,bj) = tmpFld(i,j,bi,bj) + 1.
170 ENDIF
171 ENDDO
172 kLowC(i,j,bi,bj) = 0
173 DO k= 1, Nr
174 IF ( hFacC(i,j,k,bi,bj).NE.zeroRS ) THEN
175 kLowC(i,j,bi,bj) = k
176 ENDIF
177 ENDDO
178 maskInC(i,j,bi,bj) = 0.
179 IF ( kSurfC(i,j,bi,bj).LE.Nr ) maskInC(i,j,bi,bj) = 1.
180 ENDDO
181 ENDDO
182 C- end bi,bj loops.
183 ENDDO
184 ENDDO
185
186 #ifdef ALLOW_SHELFICE
187 #ifdef ALLOW_SHELFICE_REMESHING
188 IF ( useShelfIce ) THEN
189 C-- Modify lopping factor hFacC : Remove part outside of the domain
190 C taking into account the Reference (=at rest) Surface Position Ro_shelfIce
191 CALL SHELFICE_DIG_SHELF( myThid )
192 ENDIF
193 #endif
194 #endif /* ALLOW_SHELFICE */
195
196 IF ( plotLevel.GE.debLevB ) THEN
197 c CALL PLOT_FIELD_XYRS( tmpFld,
198 c & 'Model Depths K Index' , -1, myThid )
199 CALL PLOT_FIELD_XYRS(R_low,
200 & 'Model R_low (ini_masks_etc)', -1, myThid )
201 CALL PLOT_FIELD_XYRS(Ro_surf,
202 & 'Model Ro_surf (ini_masks_etc)', -1, myThid )
203 ENDIF
204
205 C-- Calculate quantities derived from XY depth map
206 DO bj = myByLo(myThid), myByHi(myThid)
207 DO bi = myBxLo(myThid), myBxHi(myThid)
208 DO j=1-OLy,sNy+OLy
209 DO i=1-OLx,sNx+OLx
210 C Total fluid column thickness (r_unit) :
211 tmpFld(i,j,bi,bj) = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
212 C Inverse of fluid column thickness (1/r_unit)
213 IF ( tmpFld(i,j,bi,bj) .LE. zeroRS ) THEN
214 recip_Rcol(i,j,bi,bj) = 0.
215 ELSE
216 recip_Rcol(i,j,bi,bj) = 1. _d 0 / tmpFld(i,j,bi,bj)
217 ENDIF
218 ENDDO
219 ENDDO
220 ENDDO
221 ENDDO
222
223 DO bj=myByLo(myThid), myByHi(myThid)
224 DO bi=myBxLo(myThid), myBxHi(myThid)
225
226 C-- rLow & reference rSurf at Western & Southern edges (U and V points)
227 i = 1-OLx
228 DO j=1-OLy,sNy+OLy
229 rLowW (i,j,bi,bj) = rF(1)
230 rSurfW(i,j,bi,bj) = rF(1)
231 ENDDO
232 j = 1-OLy
233 DO i=1-OLx,sNx+OLx
234 rLowS (i,j,bi,bj) = rF(1)
235 rSurfS(i,j,bi,bj) = rF(1)
236 ENDDO
237 DO j=1-OLy,sNy+OLy
238 DO i=2-OLx,sNx+OLx
239 rLowW(i,j,bi,bj) =
240 & MAX( R_low(i-1,j,bi,bj), R_low(i,j,bi,bj) )
241 rSurfW(i,j,bi,bj) =
242 & MIN( Ro_surf(i-1,j,bi,bj), Ro_surf(i,j,bi,bj) )
243 rSurfW(i,j,bi,bj) =
244 & MAX( rSurfW(i,j,bi,bj), rLowW(i,j,bi,bj) )
245 ENDDO
246 ENDDO
247 DO j=2-OLy,sNy+OLy
248 DO i=1-OLx,sNx+OLx
249 rLowS(i,j,bi,bj) =
250 & MAX( R_low(i,j-1,bi,bj), R_low(i,j,bi,bj) )
251 rSurfS(i,j,bi,bj) =
252 & MIN( Ro_surf(i,j-1,bi,bj), Ro_surf(i,j,bi,bj) )
253 rSurfS(i,j,bi,bj) =
254 & MAX( rSurfS(i,j,bi,bj), rLowS(i,j,bi,bj) )
255 ENDDO
256 ENDDO
257
258 C-- hFacW and hFacS (at U and V points)
259 DO k=1, Nr
260 DO j=1-OLy,sNy+OLy
261 hFacW(1-OLx,j,k,bi,bj) = 0.
262 DO i=2-OLx,sNx+OLx
263 hFacW(i,j,k,bi,bj) =
264 & MIN( hFacC(i,j,k,bi,bj), hFacC(i-1,j,k,bi,bj) )
265 ENDDO
266 ENDDO
267 DO i=1-OLx,sNx+OLx
268 hFacS(i,1-OLy,k,bi,bj) = 0.
269 ENDDO
270 DO j=2-OLy,sNy+OLy
271 DO i=1-OLx,sNx+OLx
272 hFacS(i,j,k,bi,bj) =
273 & MIN( hFacC(i,j,k,bi,bj), hFacC(i,j-1,k,bi,bj) )
274 ENDDO
275 ENDDO
276 ENDDO
277
278 IF ( useShelfIce ) THEN
279 C Adjust reference rSurf at U and V points in order to get consistent
280 C column thickness from Sum_k(hFac*drF) and rSurf-rLow
281 DO j=1-OLy,sNy+OLy
282 DO i=1-OLx,sNx+OLx
283 rSurfW(i,j,bi,bj) = rLowW(i,j,bi,bj)
284 rSurfS(i,j,bi,bj) = rLowS(i,j,bi,bj)
285 ENDDO
286 ENDDO
287 DO k=1,Nr
288 DO j=1-OLy,sNy+OLy
289 DO i=1-OLx,sNx+OLx
290 rSurfW(i,j,bi,bj) = rSurfW(i,j,bi,bj)
291 & + hFacW(i,j,k,bi,bj)*drF(k)
292 rSurfS(i,j,bi,bj) = rSurfS(i,j,bi,bj)
293 & + hFacS(i,j,k,bi,bj)*drF(k)
294 ENDDO
295 ENDDO
296 ENDDO
297 ENDIF
298
299 C- end bi,bj loops.
300 ENDDO
301 ENDDO
302 CALL EXCH_UV_XYZ_RS( hFacW, hFacS, .FALSE., myThid )
303 CALL EXCH_UV_XY_RS( rSurfW, rSurfS, .FALSE., myThid )
304 CALL EXCH_UV_XY_RS( rLowW, rLowS, .FALSE., myThid )
305
306 C-- Addtional closing of Western and Southern grid-cell edges: for example,
307 C a) might add some "thin walls" in specific location
308 C-- b) close non-periodic N & S boundaries of lat-lon grid at the N/S poles.
309 CALL ADD_WALLS2MASKS( myThid )
310
311 C-- Calculate surface k index for interface W & S (U & V points)
312 DO bj=myByLo(myThid), myByHi(myThid)
313 DO bi=myBxLo(myThid), myBxHi(myThid)
314 DO j=1-OLy,sNy+OLy
315 DO i=1-OLx,sNx+OLx
316 kSurfW(i,j,bi,bj) = Nr+1
317 kSurfS(i,j,bi,bj) = Nr+1
318 DO k=Nr,1,-1
319 IF (hFacW(i,j,k,bi,bj).NE.zeroRS) kSurfW(i,j,bi,bj) = k
320 IF (hFacS(i,j,k,bi,bj).NE.zeroRS) kSurfS(i,j,bi,bj) = k
321 ENDDO
322 maskInW(i,j,bi,bj)= 0.
323 IF ( kSurfW(i,j,bi,bj).LE.Nr ) maskInW(i,j,bi,bj)= 1.
324 maskInS(i,j,bi,bj)= 0.
325 IF ( kSurfS(i,j,bi,bj).LE.Nr ) maskInS(i,j,bi,bj)= 1.
326 ENDDO
327 ENDDO
328 ENDDO
329 ENDDO
330
331 ELSE
332 #ifndef DISABLE_SIGMA_CODE
333 C--- Sigma and Hybrid-Sigma set-up:
334 CALL INI_SIGMA_HFAC( myThid )
335 #endif /* DISABLE_SIGMA_CODE */
336 ENDIF
337
338 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
339
340 C-- Write to disk: Total Column Thickness & hFac(C,W,S):
341 C This I/O is now done in write_grid.F
342 c CALL WRITE_FLD_XY_RS( 'Depth',' ',tmpFld,0,myThid)
343 c CALL WRITE_FLD_XYZ_RS( 'hFacC',' ',hFacC,0,myThid)
344 c CALL WRITE_FLD_XYZ_RS( 'hFacW',' ',hFacW,0,myThid)
345 c CALL WRITE_FLD_XYZ_RS( 'hFacS',' ',hFacS,0,myThid)
346
347 IF ( plotLevel.GE.debLevB ) THEN
348 CALL PLOT_FIELD_XYZRS( hFacC, 'hFacC' , Nr, 0, myThid )
349 CALL PLOT_FIELD_XYZRS( hFacW, 'hFacW' , Nr, 0, myThid )
350 CALL PLOT_FIELD_XYZRS( hFacS, 'hFacS' , Nr, 0, myThid )
351 ENDIF
352
353 C-- Masks and reciprocals of hFac[CWS]
354 DO bj = myByLo(myThid), myByHi(myThid)
355 DO bi = myBxLo(myThid), myBxHi(myThid)
356 DO k=1,Nr
357 DO j=1-OLy,sNy+OLy
358 DO i=1-OLx,sNx+OLx
359 IF ( hFacC(i,j,k,bi,bj).NE.zeroRS ) THEN
360 recip_hFacC(i,j,k,bi,bj) = 1. _d 0 / hFacC(i,j,k,bi,bj)
361 maskC(i,j,k,bi,bj) = 1.
362 ELSE
363 recip_hFacC(i,j,k,bi,bj) = 0.
364 maskC(i,j,k,bi,bj) = 0.
365 ENDIF
366 IF ( hFacW(i,j,k,bi,bj).NE.zeroRS ) THEN
367 recip_hFacW(i,j,k,bi,bj) = 1. _d 0 / hFacW(i,j,k,bi,bj)
368 maskW(i,j,k,bi,bj) = 1.
369 ELSE
370 recip_hFacW(i,j,k,bi,bj) = 0.
371 maskW(i,j,k,bi,bj) = 0.
372 ENDIF
373 IF ( hFacS(i,j,k,bi,bj).NE.zeroRS ) THEN
374 recip_hFacS(i,j,k,bi,bj) = 1. _d 0 / hFacS(i,j,k,bi,bj)
375 maskS(i,j,k,bi,bj) = 1.
376 ELSE
377 recip_hFacS(i,j,k,bi,bj) = 0.
378 maskS(i,j,k,bi,bj) = 0.
379 ENDIF
380 ENDDO
381 ENDDO
382 ENDDO
383 #ifdef NONLIN_FRSURF
384 C-- Save initial geometrical hFac factor into h0Fac (fixed in time):
385 C Note: In case 1 pkg modifies hFac (from packages_init_fixed, called
386 C later in sequence of calls) this pkg would need also to update h0Fac.
387 DO k=1,Nr
388 DO j=1-OLy,sNy+OLy
389 DO i=1-OLx,sNx+OLx
390 h0FacC(i,j,k,bi,bj) = _hFacC(i,j,k,bi,bj)
391 h0FacW(i,j,k,bi,bj) = _hFacW(i,j,k,bi,bj)
392 h0FacS(i,j,k,bi,bj) = _hFacS(i,j,k,bi,bj)
393 ENDDO
394 ENDDO
395 ENDDO
396 #endif /* NONLIN_FRSURF */
397 C- end bi,bj loops.
398 ENDDO
399 ENDDO
400
401 c #ifdef ALLOW_NONHYDROSTATIC
402 C-- Calculate "recip_hFacU" = reciprocal hfac distance/volume for W cells
403 C NOTE: not used ; computed locally in CALC_GW
404 c #endif
405
406 RETURN
407 END

  ViewVC Help
Powered by ViewVC 1.1.22