/[MITgcm]/MITgcm_contrib/shelfice_remeshing/AUTO/code/ini_masks_etc_JJ.F
ViewVC logotype

Annotation of /MITgcm_contrib/shelfice_remeshing/AUTO/code/ini_masks_etc_JJ.F

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


Revision 1.1 - (hide annotations) (download)
Thu Sep 10 14:56:34 2015 UTC (9 years, 10 months ago) by dgoldberg
Branch: MAIN
*** empty log message ***

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

  ViewVC Help
Powered by ViewVC 1.1.22