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

Annotation 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 - (hide 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 dgoldberg 1.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