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

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

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


Revision 1.2 - (hide annotations) (download)
Sat Apr 2 17:22:56 2016 UTC (9 years, 11 months ago) by dgoldberg
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +14 -1 lines
replace other digging code. split out digging code into separate S/R to be
called from ini_masks_etc

1 dgoldberg 1.2 C $Header: /u/gcmpack/MITgcm_contrib/shelfice_remeshing/DIG/code/ini_masks_etc.F,v 1.1 2016/04/01 10:19:37 dgoldberg Exp $
2 dgoldberg 1.1 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     #ifdef NONLIN_FRSURF
32     # include "SURFACE.h"
33     #endif /* NONLIN_FRSURF */
34    
35     C !INPUT/OUTPUT PARAMETERS:
36     C == Routine arguments ==
37     C myThid :: Number of this instance of INI_MASKS_ETC
38     INTEGER myThid
39    
40     C !LOCAL VARIABLES:
41     C == Local variables ==
42     C bi,bj :: tile indices
43     C i,j,k :: Loop counters
44     C tmpfld :: Temporary array used to compute & write Total Depth
45     _RS tmpfld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
46     INTEGER bi, bj
47     INTEGER i, j, k
48     _RL hFacCtmp
49     _RL hFacMnSz
50     CEOP
51    
52     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
53    
54     IF ( selectSigmaCoord.EQ.0 ) THEN
55     C--- r-coordinate with partial-cell or full cell mask
56    
57     C-- Calculate lopping factor hFacC : over-estimate the part inside of the domain
58     C taking into account the lower_R Boundary (Bathymetrie / Top of Atmos)
59     DO bj=myByLo(myThid), myByHi(myThid)
60     DO bi=myBxLo(myThid), myBxHi(myThid)
61     DO k=1, Nr
62     hFacMnSz=max( hFacMin, min(hFacMinDr*recip_drF(k),1. _d 0) )
63     DO j=1-OLy,sNy+OLy
64     DO i=1-OLx,sNx+OLx
65     C o Non-dimensional distance between grid bound. and domain lower_R bound.
66     hFacCtmp = (rF(k)-R_low(i,j,bi,bj))*recip_drF(k)
67     C o Select between, closed, open or partial (0,1,0-1)
68     hFacCtmp=min( max( hFacCtmp, 0. _d 0) , 1. _d 0)
69     C o Impose minimum fraction and/or size (dimensional)
70     IF (hFacCtmp.LT.hFacMnSz) THEN
71     IF (hFacCtmp.LT.hFacMnSz*0.5) THEN
72     hFacC(i,j,k,bi,bj)=0.
73     ELSE
74     hFacC(i,j,k,bi,bj)=hFacMnSz
75     ENDIF
76     ELSE
77     hFacC(i,j,k,bi,bj)=hFacCtmp
78     ENDIF
79     ENDDO
80     ENDDO
81     ENDDO
82    
83     C- Re-calculate lower-R Boundary position, taking into account hFacC
84     DO j=1-OLy,sNy+OLy
85     DO i=1-OLx,sNx+OLx
86     R_low(i,j,bi,bj) = rF(1)
87     ENDDO
88     ENDDO
89     DO k=Nr,1,-1
90     DO j=1-OLy,sNy+OLy
91     DO i=1-OLx,sNx+OLx
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     ENDDO
129     ENDDO
130    
131     #ifdef ALLOW_SHELFICE
132     IF ( useShelfIce ) THEN
133     C-- Modify lopping factor hFacC : Remove part outside of the domain
134     C taking into account the Reference (=at rest) Surface Position Ro_shelfIce
135     CALL SHELFICE_UPDATE_MASKS(
136     I rF, recip_drF,R_low,
137     U hFacC,
138     I myThid )
139     ENDIF
140     #endif /* ALLOW_SHELFICE */
141    
142     C- Re-calculate Reference surface position, taking into account hFacC
143     C initialize Total column fluid thickness and surface k index
144     C Note: if no fluid (continent) ==> kSurf = Nr+1
145     DO bj=myByLo(myThid), myByHi(myThid)
146     DO bi=myBxLo(myThid), myBxHi(myThid)
147     DO j=1-OLy,sNy+OLy
148     DO i=1-OLx,sNx+OLx
149     tmpfld(i,j,bi,bj) = 0.
150     kSurfC(i,j,bi,bj) = Nr+1
151     c maskH(i,j,bi,bj) = 0.
152     Ro_surf(i,j,bi,bj) = R_low(i,j,bi,bj)
153     DO k=Nr,1,-1
154     Ro_surf(i,j,bi,bj) = Ro_surf(i,j,bi,bj)
155     & + drF(k)*hFacC(i,j,k,bi,bj)
156     IF (hFacC(i,j,k,bi,bj).NE.0.) THEN
157     kSurfC(i,j,bi,bj) = k
158     c maskH(i,j,bi,bj) = 1.
159     tmpfld(i,j,bi,bj) = tmpfld(i,j,bi,bj) + 1.
160     ENDIF
161     ENDDO
162     kLowC(i,j,bi,bj) = 0
163     DO k= 1, Nr
164     IF (hFacC(i,j,k,bi,bj).NE.0) THEN
165     kLowC(i,j,bi,bj) = k
166     ENDIF
167     ENDDO
168     maskInC(i,j,bi,bj)= 0.
169     IF ( kSurfC(i,j,bi,bj).LE.Nr ) maskInC(i,j,bi,bj)= 1.
170     ENDDO
171     ENDDO
172     C- end bi,bj loops.
173     ENDDO
174     ENDDO
175    
176 dgoldberg 1.2 #ifdef ALLOW_SHELFICE
177     #ifdef SHELFICE_REMESHING
178     IF ( useShelfIce ) THEN
179     C-- Modify lopping factor hFacC : Remove part outside of the domain
180     C taking into account the Reference (=at rest) Surface Position Ro_shelfIce
181     CALL SHELFICE_DIG_SHELF(
182     U hFacC,
183     I myThid )
184     ENDIF
185     #endif
186     #endif /* ALLOW_SHELFICE */
187    
188 dgoldberg 1.1 IF ( printDomain ) THEN
189     c CALL PLOT_FIELD_XYRS( tmpfld,
190     c & 'Model Depths K Index' , -1, myThid )
191     CALL PLOT_FIELD_XYRS(R_low,
192     & 'Model R_low (ini_masks_etc)', -1, myThid )
193     CALL PLOT_FIELD_XYRS(Ro_surf,
194     & 'Model Ro_surf (ini_masks_etc)', -1, myThid )
195 dgoldberg 1.2
196 dgoldberg 1.1 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     IF ( tmpfld(i,j,bi,bj) .LE. 0. ) THEN
208     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     C-- hFacW and hFacS (at U and V points)
218     DO bj=myByLo(myThid), myByHi(myThid)
219     DO bi=myBxLo(myThid), myBxHi(myThid)
220     DO k=1, Nr
221     DO j=1-OLy,sNy+OLy
222     hFacW(1-OLx,j,k,bi,bj)= 0.
223     DO i=2-OLx,sNx+OLx
224     hFacW(i,j,k,bi,bj)=
225     & MIN(hFacC(i,j,k,bi,bj),hFacC(i-1,j,k,bi,bj))
226     ENDDO
227     ENDDO
228     DO i=1-OLx,sNx+OLx
229     hFacS(i,1-OLy,k,bi,bj)= 0.
230     ENDDO
231     DO j=2-OLy,sNy+oly
232     DO i=1-OLx,sNx+OLx
233     hFacS(i,j,k,bi,bj)=
234     & MIN(hFacC(i,j,k,bi,bj),hFacC(i,j-1,k,bi,bj))
235     ENDDO
236     ENDDO
237     ENDDO
238     C rLow & reference rSurf at Western & Southern edges (U and V points)
239     i = 1-OLx
240     DO j=1-OLy,sNy+OLy
241     rLowW (i,j,bi,bj) = 0.
242     rSurfW(i,j,bi,bj) = 0.
243     ENDDO
244     j = 1-OLy
245     DO i=1-OLx,sNx+OLx
246     rLowS (i,j,bi,bj) = 0.
247     rSurfS(i,j,bi,bj) = 0.
248     ENDDO
249     DO j=1-OLy,sNy+OLy
250     DO i=2-OLx,sNx+OLx
251     rLowW(i,j,bi,bj) =
252     & MAX( R_low(i-1,j,bi,bj), R_low(i,j,bi,bj) )
253     rSurfW(i,j,bi,bj) =
254     & MIN( Ro_surf(i-1,j,bi,bj), Ro_surf(i,j,bi,bj) )
255     rSurfW(i,j,bi,bj) =
256     & MAX( rSurfW(i,j,bi,bj), rLowW(i,j,bi,bj) )
257     ENDDO
258     ENDDO
259     DO j=2-OLy,sNy+OLy
260     DO i=1-OLx,sNx+OLx
261     rLowS(i,j,bi,bj) =
262     & MAX( R_low(i,j-1,bi,bj), R_low(i,j,bi,bj) )
263     rSurfS(i,j,bi,bj) =
264     & MIN( Ro_surf(i,j-1,bi,bj), Ro_surf(i,j,bi,bj) )
265     rSurfS(i,j,bi,bj) =
266     & MAX( rSurfS(i,j,bi,bj), rLowS(i,j,bi,bj) )
267     ENDDO
268     ENDDO
269     C- end bi,bj loops.
270     ENDDO
271     ENDDO
272     CALL EXCH_UV_XYZ_RS(hFacW,hFacS,.FALSE.,myThid)
273     CALL EXCH_UV_XY_RS( rSurfW, rSurfS, .FALSE., myThid )
274     CALL EXCH_UV_XY_RS( rLowW, rLowS, .FALSE., myThid )
275    
276     C-- Addtional closing of Western and Southern grid-cell edges: for example,
277     C a) might add some "thin walls" in specific location
278     C-- b) close non-periodic N & S boundaries of lat-lon grid at the N/S poles.
279     CALL ADD_WALLS2MASKS( myThid )
280    
281     C-- Calculate surface k index for interface W & S (U & V points)
282     DO bj=myByLo(myThid), myByHi(myThid)
283     DO bi=myBxLo(myThid), myBxHi(myThid)
284     DO j=1-OLy,sNy+OLy
285     DO i=1-OLx,sNx+OLx
286     kSurfW(i,j,bi,bj) = Nr+1
287     kSurfS(i,j,bi,bj) = Nr+1
288     DO k=Nr,1,-1
289     IF (hFacW(i,j,k,bi,bj).NE.0.) kSurfW(i,j,bi,bj) = k
290     IF (hFacS(i,j,k,bi,bj).NE.0.) kSurfS(i,j,bi,bj) = k
291     ENDDO
292     maskInW(i,j,bi,bj)= 0.
293     IF ( kSurfW(i,j,bi,bj).LE.Nr ) maskInW(i,j,bi,bj)= 1.
294     maskInS(i,j,bi,bj)= 0.
295     IF ( kSurfS(i,j,bi,bj).LE.Nr ) maskInS(i,j,bi,bj)= 1.
296     ENDDO
297     ENDDO
298     ENDDO
299     ENDDO
300    
301     ELSE
302     #ifndef DISABLE_SIGMA_CODE
303     C--- Sigma and Hybrid-Sigma set-up:
304     CALL INI_SIGMA_HFAC( myThid )
305     #endif /* DISABLE_SIGMA_CODE */
306     ENDIF
307    
308     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
309    
310     C-- Write to disk: Total Column Thickness & hFac(C,W,S):
311     C This I/O is now done in write_grid.F
312     c CALL WRITE_FLD_XY_RS( 'Depth',' ',tmpfld,0,myThid)
313     c CALL WRITE_FLD_XYZ_RS( 'hFacC',' ',hFacC,0,myThid)
314     c CALL WRITE_FLD_XYZ_RS( 'hFacW',' ',hFacW,0,myThid)
315     c CALL WRITE_FLD_XYZ_RS( 'hFacS',' ',hFacS,0,myThid)
316    
317     IF ( printDomain ) THEN
318     CALL PLOT_FIELD_XYZRS( hFacC, 'hFacC' , Nr, 0, myThid )
319     CALL PLOT_FIELD_XYZRS( hFacW, 'hFacW' , Nr, 0, myThid )
320     CALL PLOT_FIELD_XYZRS( hFacS, 'hFacS' , Nr, 0, myThid )
321     ENDIF
322    
323     C-- Masks and reciprocals of hFac[CWS]
324     DO bj = myByLo(myThid), myByHi(myThid)
325     DO bi = myBxLo(myThid), myBxHi(myThid)
326     DO k=1,Nr
327     DO j=1-OLy,sNy+OLy
328     DO i=1-OLx,sNx+OLx
329     IF (hFacC(i,j,k,bi,bj) .NE. 0. ) THEN
330     recip_hFacC(i,j,k,bi,bj) = 1. _d 0 / hFacC(i,j,k,bi,bj)
331     maskC(i,j,k,bi,bj) = 1.
332     ELSE
333     recip_hFacC(i,j,k,bi,bj) = 0.
334     maskC(i,j,k,bi,bj) = 0.
335     ENDIF
336     IF (hFacW(i,j,k,bi,bj) .NE. 0. ) THEN
337     recip_hFacW(i,j,k,bi,bj) = 1. _d 0 / hFacW(i,j,k,bi,bj)
338     maskW(i,j,k,bi,bj) = 1.
339     ELSE
340     recip_hFacW(i,j,k,bi,bj) = 0.
341     maskW(i,j,k,bi,bj) = 0.
342     ENDIF
343     IF (hFacS(i,j,k,bi,bj) .NE. 0. ) THEN
344     recip_hFacS(i,j,k,bi,bj) = 1. _d 0 / hFacS(i,j,k,bi,bj)
345     maskS(i,j,k,bi,bj) = 1.
346     ELSE
347     recip_hFacS(i,j,k,bi,bj) = 0.
348     maskS(i,j,k,bi,bj) = 0.
349     ENDIF
350     ENDDO
351     ENDDO
352     ENDDO
353     #ifdef NONLIN_FRSURF
354     C-- Save initial geometrical hFac factor into h0Fac (fixed in time):
355     C Note: In case 1 pkg modifies hFac (from packages_init_fixed, called
356     C later in sequence of calls) this pkg would need also to update h0Fac.
357     DO k=1,Nr
358     DO j=1-OLy,sNy+OLy
359     DO i=1-OLx,sNx+OLx
360     h0FacC(i,j,k,bi,bj) = _hFacC(i,j,k,bi,bj)
361     h0FacW(i,j,k,bi,bj) = _hFacW(i,j,k,bi,bj)
362     h0FacS(i,j,k,bi,bj) = _hFacS(i,j,k,bi,bj)
363     ENDDO
364     ENDDO
365     ENDDO
366     #endif /* NONLIN_FRSURF */
367     C- end bi,bj loops.
368     ENDDO
369     ENDDO
370    
371     c #ifdef ALLOW_NONHYDROSTATIC
372     C-- Calculate "recip_hFacU" = reciprocal hfac distance/volume for W cells
373     C NOTE: not used ; computed locally in CALC_GW
374     c #endif
375    
376     RETURN
377     END

  ViewVC Help
Powered by ViewVC 1.1.22