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 |