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

Annotation of /MITgcm_contrib/shelfice_remeshing/AUTO/code/shelfice_update_masks_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:35 2015 UTC (9 years, 10 months ago) by dgoldberg
Branch: MAIN
*** empty log message ***

1 dgoldberg 1.1 C $Header: /u/gcmpack/MITgcm/pkg/shelfice/shelfice_update_masks.F,v 1.6 2014/09/11 19:20:38 jmc Exp $
2     C $Name: $
3    
4     #include "SHELFICE_OPTIONS.h"
5     #ifdef ALLOW_CTRL
6     # include "CTRL_OPTIONS.h"
7     #endif
8    
9     CBOP
10     C !ROUTINE: SHELFICE_UPDATE_MASKS
11     C !INTERFACE:
12     SUBROUTINE SHELFICE_UPDATE_MASKS_JJ(
13     I rF, recip_drF,
14     U hFacC,
15     I myThid )
16     C !DESCRIPTION: \bv
17     C *==========================================================*
18     C | SUBROUTINE SHELFICE_UPDATE_MASKS
19     C | o modify topography factor hFacC according to ice shelf
20     C | topography
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 "DYNVARS.h"
31     #include "SURFACE.h"
32     #ifdef ALLOW_SHELFICE
33     # include "SHELFICE.h"
34     #endif /* ALLOW_SHELFICE */
35    
36     C !INPUT/OUTPUT PARAMETERS:
37     C == Routine arguments ==
38     C rF :: R-coordinate of face of cell (units of r).
39     C recip_drF :: Recipricol of cell face separation along Z axis ( units of r ).
40     C hFacC :: Fraction of cell in vertical which is open (see GRID.h)
41     C myThid :: Number of this instance of SHELFICE_UPDATE_MASKS
42     _RS rF (1:Nr+1)
43     _RS recip_drF (1:Nr)
44     _RS hFacC (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:Nr,nSx,nSy)
45    
46     INTEGER myThid
47    
48     #ifdef ALLOW_SHELFICE
49     C !LOCAL VARIABLES:
50     C == Local variables ==
51     C bi,bj :: tile indices
52     C I,J,K :: Loop counters
53     INTEGER bi, bj
54     INTEGER I, J, K
55     _RL hFacCtmp
56     _RL hFacMnSz
57     c _RS R_shelfIceOld (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
58     CEOP
59    
60    
61     C initialize R_shelfIce
62     c DO bj = myByLo(myThid), myByHi(myThid)
63     c DO bi = myBxLo(myThid), myBxHi(myThid)
64     c DO j=1-OLy,sNy+OLy
65     c DO i=1-OLx,sNx+OLx
66     c R_shelfIceOld(i,j,bi,bj) = R_shelfIce(i,j,bi,bj)
67     c R_shelfIce(i,j,bi,bj) = 0. _d 0
68     c ENDDO
69     c ENDDO
70     c ENDDO
71     c ENDDO
72    
73    
74     c IF ( SHELFICEtopoFile .NE. ' ' ) THEN
75     c _BARRIER
76     C Read the shelfIce draught using the mid-level I/O pacakage read_write_rec
77     C The 0 is the "iteration" argument. The 1 is the record number.
78     c CALL READ_REC_XY_RS( SHELFICEtopoFile, R_shelfIce,
79     c & 1, 0, myThid )
80     C Read the shelfIce draught using the mid-level I/O pacakage read_write_fld
81     C The 0 is the "iteration" argument. The ' ' is an empty suffix
82     C CALL READ_FLD_XY_RS( SHELFICEtopoFile, ' ', R_shelfIce,
83     C & 0, myThid )
84     c ELSE
85     c DO bj = myByLo(myThid), myByHi(myThid)
86     c DO bi = myBxLo(myThid), myBxHi(myThid)
87     c DO j=1-OLy,sNy+OLy
88     c DO i=1-OLx,sNx+OLx
89     c
90     c R_shelfIce(i,j,bi,bj) = R_shelfIceOld(i,j,bi,bj)
91     c
92     c ENDDO
93     c ENDDO
94     c ENDDO
95     c ENDDO
96    
97     c ENDIF
98    
99    
100     C- Update etaN
101     DO bj = myByLo(myThid), myByHi(myThid)
102     DO bi = myBxLo(myThid), myBxHi(myThid)
103     DO J = 1-OLy,sNy+OLy
104     DO I = 1-OLx,sNx+OLx
105     IF ( R_shelfice(I,J,bi,bj) .LT. 0.0) THEN
106     IF (etah(I,J,bi,bj) .GT. SHELFICESplitThreshold ) THEN
107     K = MAX(1,kTopC(I,J,bi,bj))
108     etaN(I,J,bi,bj) = etaN(I,J,bi,bj) - 1/recip_drF(K)
109     etaH(I,J,bi,bj) = etaH(I,J,bi,bj) - 1/recip_drF(K)
110     R_shelfIce(I,J,bi,bj) = R_shelfIce(I,J,bi,bj)+1/recip_drF(K)
111     uVel(I,J,K-1,bi,bj)=uVel(I,J,K,bi,bj)
112     uVel(I+1,J,K-1,bi,bj)=uVel(I+1,J,K,bi,bj)
113     vVel(I,J,K-1,bi,bj)=vVel(I,J,K,bi,bj)
114     vVel(I,J+1,K-1,bi,bj)=vVel(I,J+1,K,bi,bj)
115     gvnm1(I,J,K-1,bi,bj)=0.0
116     gvnm1(I,J+1,K-1,bi,bj)=0.0
117     gunm1(I,J,K-1,bi,bj)=0.0
118     gunm1(I+1,J,K-1,bi,bj)=0.0
119     salt(I,J,K-1,bi,bj)=salt(I,J,K,bi,bj)
120     theta(I,J,K-1,bi,bj)=theta(I,J,K,bi,bj)
121    
122     hfac_surfc(I,J,bi,bj)=hfac_surfc(I,J,bi,bj)-1.0
123     hfac_surfnm1c(I,J,bi,bj)=hfac_surfc(I,J,bi,bj)
124     hfacC(I,J,K,bi,bj)=1.0
125    
126     ENDIF
127     IF (etah(I,J,bi,bj) .LT. SHELFICEMergeThreshold ) THEN
128     K = MAX(1,kTopC(I,J,bi,bj))
129    
130     salt(I,J,K+1,bi,bj)=((salt(I,J,K,bi,bj)*(1/recip_drF(K)+
131     & etaN(I,J,bi,bj)))+(salt(I,J,K+1,bi,bj)*1/recip_drF(K+1)))/(
132     & 1/recip_drF(K)+1/recip_drF(K+1)+etaN(I,J,bi,bj))
133    
134     theta(I,J,K+1,bi,bj)=((theta(I,J,K,bi,bj)*(1/recip_drF(K)+
135     & etaN(I,J,bi,bj)))+(theta(I,J,K+1,bi,bj)*1/recip_drF(K+1)))/(
136     & 1/recip_drF(K)+1/recip_drF(K+1)+etaN(I,J,bi,bj))
137    
138     vVel(I,J,K+1,bi,bj)=((vVel(I,J,K,bi,bj)*(1/recip_drF(K)+
139     & etaN(I,J,bi,bj)))+(vVel(I,J,K+1,bi,bj)*1/recip_drF(K+1)))/(
140     & 1/recip_drF(K)+1/recip_drF(K+1)+etaN(I,J,bi,bj))
141    
142     vVel(I,J+1,K+1,bi,bj)=((vVel(I,J+1,K,bi,bj)*(1/recip_drF(K)+
143     & etaN(I,J,bi,bj)))+(vVel(I,J+1,K+1,bi,bj)*1/recip_drF(K+1)))/
144     & (1/recip_drF(K)+1/recip_drF(K+1)+etaN(I,J,bi,bj))
145    
146     uVel(I,J,K+1,bi,bj)=((uVel(I,J,K,bi,bj)*(1/recip_drF(K)+
147     & etaN(I,J,bi,bj)))+(uVel(I,J,K+1,bi,bj)*1/recip_drF(K+1)))/(
148     & 1/recip_ drF(K)+1/recip_drF(K+1)+etaN(I,J,bi,bj))
149    
150     uVel(I+1,J,K+1,bi,bj)=((uVel(I+1,J,K,bi,bj)*(1/recip_drF(K)+
151     & etaN(I,J,bi,bj)))+(uVel(I+1,J,K+1,bi,bj)*1/recip_drF(K+1)))/
152     & (1/recip_drF(K)+1/recip_drF(K+1)+etaN(I,J,bi,bj))
153    
154     etaN(I,J,bi,bj) = etaN(I,J,bi,bj) +1/recip_drF(K)
155     etaH(I,J,bi,bj) = etaH(I,J,bi,bj) +1/recip_drF(K)
156     R_shelfice(I,J,bi,bj) = R_shelfice(I,J,bi,bj) -1/recip_drF(K)
157    
158     gvnm1(I,J,K+1,bi,bj)=0.0
159     gvnm1(I,J+1,K+1,bi,bj)=0.0
160     gunm1(I,J,K+1,bi,bj)=0.0
161     gunm1(I+1,J,K+1,bi,bj)=0.0
162    
163     hfac_surfc(I,J,bi,bj)=hfac_surfc(I,J,bi,bj)+1.0
164     c hfac_surfnm1c(I,J,bi,bj)=hfac_surfc(I,J,bi,bj)
165     hfacC(I,J,K,bi,bj)=1.0
166     ENDIF
167     ENDIF
168     ENDDO
169     ENDDO
170     ENDDO
171     ENDDO
172    
173     DO bj = myByLo(myThid), myByHi(myThid)
174     DO bi = myBxLo(myThid), myBxHi(myThid)
175     DO J = 1-OLy,sNy+OLy
176     DO I = 1-OLx,sNx+OLx
177     etaH(I,J,bi,bj)=etaN(I,J,bi,bj)
178     etaHnm1(I,J,bi,bj)=etaH(I,J,bi,bj)
179     ENDDO
180     ENDDO
181     ENDDO
182     ENDDO
183    
184     DO bj = myByLo(myThid), myByHi(myThid)
185     DO bi = myBxLo(myThid), myBxHi(myThid)
186     DO J = 1-OLy,sNy+OLy
187     DO I = 1-OLx,sNx+OLx
188     DO K = 1, Nr
189     c gvnm1(I,J,K,bi,bj)=0.0
190     c gunm1(I,J,K,bi,bj)=0.0
191     c gwnm1(I,J,K,bi,bj)=0.0
192     ENDDO
193     ENDDO
194     ENDDO
195     ENDDO
196     ENDDO
197    
198    
199    
200     C- fill in the overlap (+ BARRIER):
201     _EXCH_XY_RS(R_shelfIce, myThid )
202    
203     C-- Calculate lopping factor hFacC : Remove part outside of the domain
204     C taking into account the Reference (=at rest) Surface Position Ro_shelfIce
205     DO bj=myByLo(myThid), myByHi(myThid)
206     DO bi=myBxLo(myThid), myBxHi(myThid)
207    
208     C-- compute contributions of shelf ice to looping factors
209     DO K=1, Nr
210     hFacMnSz=max( hFacMin, min(hFacMinDr*recip_drF(k),1. _d 0) )
211     DO J=1-OLy,sNy+OLy
212     DO I=1-OLx,sNx+OLx
213     C o Non-dimensional distance between grid boundary and model surface
214     hFacCtmp = (rF(k)-R_shelfIce(I,J,bi,bj))*recip_drF(K)
215     C o Reduce the previous fraction : substract the outside part.
216     hFacCtmp = hFacC(I,J,K,bi,bj) - max( hFacCtmp, 0. _d 0)
217     C o set to zero if empty Column :
218     hFacCtmp = max( hFacCtmp, 0. _d 0)
219     C o Impose minimum fraction and/or size (dimensional)
220     IF (hFacCtmp.LT.hFacMnSz) THEN
221     IF (hFacCtmp.LT.hFacMnSz*0.5) THEN
222     hFacC(I,J,K,bi,bj)=0.
223     ELSE
224     hFacC(I,J,K,bi,bj)=hFacMnSz
225     ENDIF
226     ELSE
227     hFacC(I,J,K,bi,bj)=hFacCtmp
228     ENDIF
229     ENDDO
230     ENDDO
231     ENDDO
232    
233     #ifdef ALLOW_SHIFWFLX_CONTROL
234     C maskSHI is a hack to play along with the general ctrl-package
235     C infrastructure, where only the k=1 layer of a 3D mask is used
236     C for 2D fields. We cannot use maskInC instead, because routines
237     C like ctrl_get_gen and ctrl_set_unpack_xy require 3D masks.
238     DO K=1,Nr
239     DO J=1-OLy,sNy+OLy
240     DO I=1-OLx,sNx+OLx
241     maskSHI(I,J,K,bi,bj) = 0. _d 0
242     ENDDO
243     ENDDO
244     ENDDO
245     DO K=1,Nr
246     DO J=1-OLy,sNy+OLy
247     DO I=1-OLx,sNx+OLx
248     IF ( ABS(R_shelfice(I,J,bi,bj)) .GT. 0. _d 0
249     & .AND. hFacC(I,J,K,bi,bj) .NE. 0. _d 0 ) THEN
250     maskSHI(I,J,K,bi,bj) = 1. _d 0
251     maskSHI(I,J,1,bi,bj) = 1. _d 0
252     ENDIF
253     ENDDO
254     ENDDO
255     ENDDO
256     #endif /* ALLOW_SHIFWFLX_CONTROL */
257    
258     C - end bi,bj loops.
259     ENDDO
260     ENDDO
261     #endif /* ALLOW_SHELFICE */
262     RETURN
263     END

  ViewVC Help
Powered by ViewVC 1.1.22