/[MITgcm]/MITgcm_contrib/verification_other/shelfice_remeshing/code/ini_masks_etc.F
ViewVC logotype

Annotation of /MITgcm_contrib/verification_other/shelfice_remeshing/code/ini_masks_etc.F

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


Revision 1.2 - (hide annotations) (download)
Tue Apr 4 23:34:37 2017 UTC (8 years, 3 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint66f
Changes since 1.1: +78 -57 lines
bring up to date with standard version in model/src

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

  ViewVC Help
Powered by ViewVC 1.1.22