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

Annotation of /MITgcm_contrib/shelfice_remeshing/CLEAN/code/shelfice_update_masks_JJ.F

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


Revision 1.5 - (hide annotations) (download)
Fri Jan 22 10:26:50 2016 UTC (9 years, 6 months ago) by dgoldberg
Branch: MAIN
CVS Tags: HEAD
Changes since 1.4: +21 -12 lines
Overlap bug fixing

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

  ViewVC Help
Powered by ViewVC 1.1.22