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

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

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


Revision 1.2 - (hide annotations) (download)
Mon Sep 28 12:59:59 2015 UTC (9 years, 9 months ago) by dgoldberg
Branch: MAIN
Changes since 1.1: +65 -22 lines
*** empty log message ***

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

  ViewVC Help
Powered by ViewVC 1.1.22