/[MITgcm]/MITgcm_contrib/dgoldberg/CPL1/code/shelfice_update_masks_remesh.F
ViewVC logotype

Annotation of /MITgcm_contrib/dgoldberg/CPL1/code/shelfice_update_masks_remesh.F

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


Revision 1.1 - (hide annotations) (download)
Wed Jul 6 18:01:26 2016 UTC (9 years ago) by dgoldberg
Branch: MAIN
CVS Tags: HEAD
moving experiment out of shelfice_remeshing to replace with vertical remeshing only

1 dgoldberg 1.1 C $Header: /u/gcmpack/MITgcm_contrib/verification_other/shelfice_remeshing/code/shelfice_update_masks_remesh.F,v 1.2 2016/05/25 13:10:42 dgoldberg 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_REMESH(
13     I rF, recip_drF, drF, kLowC,
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 drF (1:Nr)
45     _RS hFacC (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:Nr,nSx,nSy)
46     INTEGER kLowC (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
47    
48     INTEGER myThid
49    
50     #ifdef ALLOW_SHELFICE
51     #ifdef ALLOW_SHELFICE_REMESHING
52     C !LOCAL VARIABLES:
53     C == Local variables ==
54     C bi,bj :: tile indices
55     C I,J,K :: Loop counters
56     INTEGER bi, bj
57     INTEGER I, J, K
58     _RL hFacCtmp
59     _RL hFacMnSz
60    
61    
62    
63    
64     C- Update etaN
65     DO bj = myByLo(myThid), myByHi(myThid)
66     DO bi = myBxLo(myThid), myBxHi(myThid)
67     DO J = 1,sNy
68     DO I = 1,sNx
69     IF ( R_shelfice(I,J,bi,bj) .LT. 0.0) THEN
70     K = MAX(1,kTopC(I,J,bi,bj))
71     IF (etah(I,J,bi,bj)*recip_drF(K)+1.0
72     & .GT. SHELFICESplitThreshold ) THEN
73     IF (etah(I,J,bi,bj)*recip_drF(K+1)+1.0
74     & .GT. SHELFICESplitThreshold ) THEN
75     etaN(I,J,bi,bj) = etaN(I,J,bi,bj)- drF(K-1)
76     etaH(I,J,bi,bj) = etaH(I,J,bi,bj)- drF(K-1)
77     R_shelfIce(I,J,bi,bj) = R_shelfIce(I,J,bi,bj)+drF(K-1)
78     ! Rmin_surf(I,J,bi,bj) = Rmin_surf(I,J,bi,bj)+
79     ! & drF(K-1)*hFacC(I,J,K-1,bi,bj)
80    
81     uVel(I,J,K-1,bi,bj)=uVel(I,J,K,bi,bj)
82     uVel(I+1,J,K-1,bi,bj)=uVel(I+1,J,K,bi,bj)
83     vVel(I,J,K-1,bi,bj)=vVel(I,J,K,bi,bj)
84     vVel(I,J+1,K-1,bi,bj)=vVel(I,J+1,K,bi,bj)
85     gvnm1(I,J,K-1,bi,bj)=gvnm1(I,J,K,bi,bj)
86     gvnm1(I,J+1,K-1,bi,bj)=gvnm1(I,J+1,K,bi,bj)
87     gunm1(I,J,K-1,bi,bj)=gunm1(I,J,K,bi,bj)
88     gunm1(I+1,J,K-1,bi,bj)=gunm1(I,J,K,bi,bj)
89    
90     salt(I,J,K-1,bi,bj)=salt(I,J,K,bi,bj)
91     theta(I,J,K-1,bi,bj)=theta(I,J,K,bi,bj)
92     ENDIF
93     ENDIF
94     IF (kTopC(i,j,bi,bj) .LT.kLowC (i,j,bi,bj))THEN
95     K = MAX(1,kTopC(I,J,bi,bj))
96     IF (etah(I,J,bi,bj)*recip_drF(K)+1.0 .LT.
97     & SHELFICEMergeThreshold ) THEN
98     IF (etah(I,J,bi,bj)*recip_drF(K+1)+1.0 .LT.
99     & SHELFICEMergeThreshold ) THEN
100    
101     etaN(I,J,bi,bj) = etaN(I,J,bi,bj) +drF(K+1)*
102     & hFacC(i,j,K+1,bi,bj)
103     etaH(I,J,bi,bj) = etaH(I,J,bi,bj) +drF(K+1)*
104     & hFacC(i,j,K+1,bi,bj)
105     R_shelfice(I,J,bi,bj) = R_shelfice(I,J,bi,bj) -drF(K+1)*
106     & hFacC(i,j,K+1,bi,bj)
107     ! Rmin_surf(I,J,bi,bj) = Rmin_surf(I,J,bi,bj) -drF(K+1)*
108     ! & hFacC(i,j,K+1,bi,bj)
109    
110     vVel(I,J,K+1,bi,bj)=((vVel(I,J,K,bi,bj)*(drF(K)+
111     & etaN(I,J,bi,bj)))+(vVel(I,J,K+1,bi,bj)*drF(K+1)*
112     & hFacC(i,j,K+1,bi,bj)))/(drF(K)+drF(K+1)*
113     & hFacC(i,j,K+1,bi,bj)+etaN(I,J,bi,bj))
114    
115     vVel(I,J+1,K+1,bi,bj)=((vVel(I,J+1,K,bi,bj)*(drF(K)+
116     & etaN(I,J,bi,bj)))+(vVel(I,J+1,K+1,bi,bj)*drF(K+1)*
117     & hFacC(i,j,K+1,bi,bj)))/(drF(K)+drF(K+1)*
118     & hFacC(i,j,K+1,bi,bj)+etaN(I,J,bi,bj))
119    
120     uVel(I,J,K+1,bi,bj)=((uVel(I,J,K,bi,bj)*(drF(K)+
121     & etaN(I,J,bi,bj)))+(uVel(I,J,K+1,bi,bj)*drF(K+1)*
122     & hFacC(i,j,K+1,bi,bj)))/(drF(K)+drF(K+1)*
123     & hFacC(i,j,K+1,bi,bj)+etaN(I,J,bi,bj))
124    
125     uVel(I+1,J,K+1,bi,bj)=((uVel(I+1,J,K,bi,bj)*(drF(K)+
126     & etaN(I,J,bi,bj)))+(uVel(I+1,J,K+1,bi,bj)*drF(K+1)*
127     & hFacC(i,j,K+1,bi,bj)))/(drF(K)+drF(K+1)*
128     & hFacC(i,j,K+1,bi,bj)+etaN(I,J,bi,bj))
129    
130     gunm1(I+1,J,K+1,bi,bj)=((gunm1(I+1,J,K,bi,bj)*(drF(K)+
131     & etaN(I,J,bi,bj)))+(gunm1(I+1,J,K+1,bi,bj)*drF(K+1)*
132     & hFacC(i,j,K+1,bi,bj)))/(drF(K)+drF(K+1)*
133     & hFacC(i,j,K+1,bi,bj)+etaN(I,J,bi,bj))
134    
135     gunm1(I,J,K+1,bi,bj)=((gunm1(I,J,K,bi,bj)*(drF(K)+
136     & etaN(I,J,bi,bj)))+(gunm1(I,J,K+1,bi,bj)*drF(K+1)*
137     & hFacC(i,j,K+1,bi,bj)))/(drF(K)+drF(K+1)*
138     & hFacC(i,j,K+1,bi,bj)+etaN(I,J,bi,bj))
139    
140     gvnm1(I,J+1,K+1,bi,bj)=((gvnm1(I,J+1,K,bi,bj)*(drF(K)+
141     & etaN(I,J,bi,bj)))+(gvnm1(I,J+1,K+1,bi,bj)*drF(K+1)*
142     & hFacC(i,j,K+1,bi,bj)))/(drF(K)+drF(K+1)*
143     & hFacC(i,j,K+1,bi,bj)+etaN(I,J,bi,bj))
144    
145     gvnm1(I,J,K+1,bi,bj)=((gvnm1(I,J,K,bi,bj)*(drF(K)+
146     & etaN(I,J,bi,bj)))+(gvnm1(I,J,K+1,bi,bj)*drF(K+1)*
147     & hFacC(i,j,K+1,bi,bj)))/(drF(K)+drF(K+1)*
148     & hFacC(i,j,K+1,bi,bj)+etaN(I,J,bi,bj))
149    
150     salt(I,J,K+1,bi,bj)=((salt(I,J,K,bi,bj)*(drF(K)+
151     & etaN(I,J,bi,bj)))+(salt(I,J,K+1,bi,bj)*drF(K+1)*
152     & hFacC(i,j,K+1,bi,bj)))/(drF(K)+drF(K+1)*
153     & hFacC(i,j,K+1,bi,bj)+etaN(I,J,bi,bj))
154    
155     theta(I,J,K+1,bi,bj)=((theta(I,J,K,bi,bj)*(drF(K)+
156     & etaN(I,J,bi,bj)))+(theta(I,J,K+1,bi,bj)*drF(K+1)*
157     & hFacC(i,j,K+1,bi,bj)))/(drF(K)+drF(K+1)*
158     & hFacC(i,j,K+1,bi,bj)+etaN(I,J,bi,bj))
159     ENDIF
160     ENDIF
161     ENDIF
162     ENDIF
163     ENDDO
164     ENDDO
165     ENDDO
166     ENDDO
167    
168     DO bj = myByLo(myThid), myByHi(myThid)
169     DO bi = myBxLo(myThid), myBxHi(myThid)
170     DO J = 1,sNy
171     DO I = 1,sNx
172     etaH(I,J,bi,bj)=etaN(I,J,bi,bj)
173     etaHnm1(I,J,bi,bj)=etaH(I,J,bi,bj)
174     ENDDO
175     ENDDO
176     ENDDO
177     ENDDO
178    
179     DO bj = myByLo(myThid), myByHi(myThid)
180     DO bi = myBxLo(myThid), myBxHi(myThid)
181     DO J = 1,sNy
182     DO I = 1,sNx
183     K = MAX(1,kTopC(I,J,bi,bj))
184     hfac_surfc(I,J,bi,bj)= ((etaH(I,J,bi,bJ) +(drF(K)))
185     & *recip_drF(K))
186     ENDDO
187     ENDDO
188     ENDDO
189     ENDDO
190    
191     CALL EXCH_XYZ_RL(salt,myThid)
192     CALL EXCH_XYZ_RL(theta,myThid)
193     CALL EXCH_XYZ_RL(uVel,myThid)
194     CALL EXCH_XYZ_RL(vVel,myThid)
195     CALL EXCH_XYZ_RL(gunm1,myThid)
196     CALL EXCH_XYZ_RL(gvnm1,myThid)
197     CALL EXCH_XYZ_RL(hFacC,myThid)
198    
199     CALL EXCH_XY_RL(EtaN,myThid)
200     CALL EXCH_XY_RL(EtaH,myThid)
201     CALL EXCH_XY_RL(EtaHnm1,myThid)
202     CALL EXCH_XY_RL(R_shelfice,myThid)
203     ! CALL EXCH_XY_RL(Rmin_surf,myThid)
204     CALL EXCH_XY_RL(hFac_surfC,myThid)
205    
206    
207     C- fill in the overlap (+ BARRIER):
208     _EXCH_XY_RS(R_shelfIce, myThid )
209    
210     C-- Calculate lopping factor hFacC : Remove part outside of the domain
211     C taking into account the Reference (=at rest) Surface Position Ro_shelfIce
212     DO bj=myByLo(myThid), myByHi(myThid)
213     DO bi=myBxLo(myThid), myBxHi(myThid)
214    
215     C-- compute contributions of shelf ice to looping factors
216     DO K=1, Nr
217     hFacMnSz=max( hFacMin, min(hFacMinDr*recip_drF(k),1. _d 0) )
218     DO J=1-OLy,sNy+OLy
219     DO I=1-OLx,sNx+OLx
220     C o Non-dimensional distance between grid boundary and model surface
221     hFacCtmp = (rF(k)-R_shelfIce(I,J,bi,bj))*recip_drF(K)
222     C o Reduce the previous fraction : substract the outside part.
223     hFacCtmp = hFacC(I,J,K,bi,bj) - max( hFacCtmp, 0. _d 0)
224     C o set to zero if empty Column :
225     hFacCtmp = max( hFacCtmp, 0. _d 0)
226     C o Impose minimum fraction and/or size (dimensional)
227     IF (hFacCtmp.LT.hFacMnSz) THEN
228     IF (hFacCtmp.LT.hFacMnSz*0.5) THEN
229     hFacC(I,J,K,bi,bj)=0.
230     ELSE
231     hFacC(I,J,K,bi,bj)=hFacMnSz
232     ENDIF
233     ELSE
234     hFacC(I,J,K,bi,bj)=hFacCtmp
235     ENDIF
236     ENDDO
237     ENDDO
238     ENDDO
239    
240     #ifdef ALLOW_SHIFWFLX_CONTROL
241     C maskSHI is a hack to play along with the general ctrl-package
242     C infrastructure, where only the k=1 layer of a 3D mask is used
243     C for 2D fields. We cannot use maskInC instead, because routines
244     C like ctrl_get_gen and ctrl_set_unpack_xy require 3D masks.
245     DO K=1,Nr
246     DO J=1-OLy,sNy+OLy
247     DO I=1-OLx,sNx+OLx
248     maskSHI(I,J,K,bi,bj) = 0. _d 0
249     ENDDO
250     ENDDO
251     ENDDO
252     DO K=1,Nr
253     DO J=1-OLy,sNy+OLy
254     DO I=1-OLx,sNx+OLx
255     IF ( ABS(R_shelfice(I,J,bi,bj)) .GT. 0. _d 0
256     & .AND. hFacC(I,J,K,bi,bj) .NE. 0. _d 0 ) THEN
257     maskSHI(I,J,K,bi,bj) = 1. _d 0
258     maskSHI(I,J,1,bi,bj) = 1. _d 0
259     ENDIF
260     ENDDO
261     ENDDO
262     ENDDO
263     #endif /* ALLOW_SHIFWFLX_CONTROL */
264    
265     C - end bi,bj loops.
266     ENDDO
267     ENDDO
268     #endif /* ALLOW_SHELFICE */
269     #endif
270     RETURN
271     END

  ViewVC Help
Powered by ViewVC 1.1.22