1 |
C $Header: /u/gcmpack/MITgcm/pkg/streamice/streamice_adv_flux_fl_x.F,v 1.2 2013/07/02 19:09:42 jmc Exp $ |
2 |
C $Name: $ |
3 |
|
4 |
#include "STREAMICE_OPTIONS.h" |
5 |
|
6 |
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
7 |
|
8 |
SUBROUTINE STREAMICE_ADV_FLUX_FL_X ( myThid , |
9 |
I UADV , |
10 |
I TRAC , |
11 |
I BC_FACEMASK, |
12 |
I BC_XVALUES, |
13 |
O XFLUX, |
14 |
I time_step ) |
15 |
|
16 |
IMPLICIT NONE |
17 |
|
18 |
C O hflux_x ! flux per unit width across face |
19 |
C O h |
20 |
C I time_step |
21 |
|
22 |
C === Global variables === |
23 |
#include "SIZE.h" |
24 |
#include "GRID.h" |
25 |
#include "EEPARAMS.h" |
26 |
#include "PARAMS.h" |
27 |
#include "STREAMICE.h" |
28 |
!#include "GAD_FLUX_LIMITER.h" |
29 |
|
30 |
INTEGER myThid |
31 |
_RL UADV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
32 |
_RL TRAC (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
33 |
_RS BC_FACEMASK (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
34 |
_RL BC_XVALUES (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
35 |
_RL XFLUX (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
36 |
_RL time_step |
37 |
|
38 |
#ifdef ALLOW_STREAMICE |
39 |
|
40 |
C LOCAL VARIABLES |
41 |
|
42 |
INTEGER i, j, bi, bj, Gi, Gj, k |
43 |
_RL uface, phi, cfl, Cr, rdenom, d0, d1, psi |
44 |
_RL stencil (-1:1) |
45 |
LOGICAL H0_valid(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
46 |
! there are valid cells to calculate a |
47 |
! slope-limited 2nd order flux |
48 |
_RL SLOPE_LIMITER |
49 |
! _RL total_vol_out |
50 |
external SLOPE_LIMITER |
51 |
|
52 |
! total_vol_out = 0.0 |
53 |
|
54 |
DO bj=myByLo(myThid),myByHi(myThid) |
55 |
DO bi=myBxLo(myThid),myBxHi(myThid) |
56 |
DO j=1-oly,sNy+oly |
57 |
DO i=1-olx,sNx+olx |
58 |
H0_valid(i,j,bi,bj)=.false. |
59 |
ENDDO |
60 |
ENDDO |
61 |
ENDDO |
62 |
ENDDO |
63 |
|
64 |
DO bj=myByLo(myThid),myByHi(myThid) |
65 |
DO bi=myBxLo(myThid),myBxHi(myThid) |
66 |
DO j=1-3,sNy+3 |
67 |
Gj = (myYGlobalLo-1)+(bj-1)*sNy+j |
68 |
IF (((Gj .ge. 1) .and. (Gj .le. Ny)) |
69 |
& .or.STREAMICE_NS_PERIODIC) THEN |
70 |
DO i=1,sNx+1 |
71 |
C THESE ARRAY BOUNDS INSURE THAT AFTER THIS STEP, |
72 |
C VALUES WILL BE RELIABLE 1 GRID CELLS OUT IN THE |
73 |
C X DIRECTION AND 3 CELLS OUT IN THE Y DIR |
74 |
IF ((STREAMICE_hmask(i,j,bi,bj).eq.1.0) .or. |
75 |
& ((STREAMICE_hmask(i-1,j,bi,bj).eq.1.0) .and. |
76 |
& (STREAMICE_hmask(i,j,bi,bj).ne.1.0))) THEN |
77 |
|
78 |
Gi = (myXGlobalLo-1)+(bi-1)*sNx+i |
79 |
|
80 |
uface = UADV(i,j,bi,bj) |
81 |
cfl = ABS(uface) * time_step * recip_dxC(i,j,bi,bj) |
82 |
|
83 |
IF (BC_FACEMASK(i,j,bi,bj).eq.3.0 .and. |
84 |
& uface.gt.0 .and. |
85 |
& STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN |
86 |
XFLUX (i,j,bi,bj) = BC_XVALUES(i,j,bi,bj) * uface |
87 |
ELSEIF |
88 |
& (BC_FACEMASK(i,j,bi,bj).eq.3.0 .and. |
89 |
& uface.le.0 .and. |
90 |
& STREAMICE_hmask(i-1,j,bi,bj).eq.1.0) THEN |
91 |
XFLUX (i,j,bi,bj) = BC_XVALUES(i,j,bi,bj) * uface |
92 |
ELSE |
93 |
|
94 |
IF (uface .gt. 0. _d 0) THEN |
95 |
DO k=-1,1 |
96 |
stencil (k) = TRAC(i+k-1,j,bi,bj) |
97 |
ENDDO |
98 |
IF ((STREAMICE_hmask(i,j,bi,bj).eq.1.0) .and. |
99 |
& (STREAMICE_hmask(i-2,j,bi,bj).eq.1.0)) |
100 |
& H0_valid(i,j,bi,bj)=.true. |
101 |
|
102 |
IF (((Gi.eq.1).and.(STREAMICE_hmask(i-1,j,bi,bj).eq.3.0)) |
103 |
& .and.(.not.STREAMICE_EW_PERIODIC)) |
104 |
& THEN ! we are at western bdry and there is a thick. bdry cond |
105 |
|
106 |
XFLUX (i,j,bi,bj) = TRAC(i-1,j,bi,bj) * uface |
107 |
|
108 |
ELSEIF (H0_valid(i,j,bi,bj)) THEN |
109 |
|
110 |
rdenom = (stencil(1)-stencil(0)) |
111 |
IF (rdenom .ne. 0.) THEN |
112 |
Cr = (stencil(0)-stencil(-1))/rdenom |
113 |
ELSE |
114 |
Cr = 1.E20 * (stencil(0)-stencil(-1)) |
115 |
ENDIF |
116 |
|
117 |
IF (STREAMICE_ADV_SCHEME.ne.'DST3') THEN |
118 |
! phi = SLOPE_LIMITER(stencil(0)-stencil(-1), |
119 |
! & stencil(1)-stencil(0)) |
120 |
phi = SLOPE_LIMITER (Cr) |
121 |
ELSE |
122 |
d0 = (2.-cfl)*(1.-cfl)/6.0 |
123 |
d1 = (1.-cfl**2)/6.0 |
124 |
psi = d0+d1*Cr |
125 |
phi = MAX(0. _d 0,MIN(MIN(1. _d 0,psi), |
126 |
& Cr*(1. _d 0 -CFL)/(CFL+1. _d -20) )) |
127 |
ENDIF |
128 |
|
129 |
IF (STREAMICE_ADV_SCHEME.ne.'DST3') THEN |
130 |
XFLUX (i,j,bi,bj) = uface * |
131 |
& (stencil(0) + phi * .5 * (1.0-cfl) * |
132 |
& (stencil(1)-stencil(0))) |
133 |
ELSE |
134 |
XFLUX (i,j,bi,bj) = uface * |
135 |
& (stencil(0) + phi * |
136 |
& (stencil(1)-stencil(0))) |
137 |
ENDIF |
138 |
|
139 |
ELSE ! one of the two cells needed for a HO scheme is missing, use FO scheme |
140 |
|
141 |
XFLUX (i,j,bi,bj) = uface * stencil(0) |
142 |
|
143 |
ENDIF |
144 |
|
145 |
ELSEIF (uface .lt. 0. _d 0) THEN ! uface <= 0 |
146 |
|
147 |
DO k=-1,1 |
148 |
stencil (k) = TRAC(i-k,j,bi,bj) |
149 |
ENDDO |
150 |
IF ((STREAMICE_hmask(i-1,j,bi,bj).eq.1.0) .and. |
151 |
& (STREAMICE_hmask(i+1,j,bi,bj).eq.1.0)) |
152 |
& H0_valid(i,j,bi,bj)=.true. |
153 |
|
154 |
IF (((Gi.eq.Nx).and.(STREAMICE_hmask(i+1,j,bi,bj).eq.3.0)) |
155 |
& .and.(.not.STREAMICE_EW_PERIODIC)) |
156 |
& THEN ! we are at western bdry and there is a thick. bdry cond |
157 |
|
158 |
XFLUX (i,j,bi,bj) = TRAC(i+1,j,bi,bj) * uface |
159 |
|
160 |
ELSEIF (H0_valid(i,j,bi,bj)) THEN |
161 |
|
162 |
rdenom = (stencil(1)-stencil(0)) |
163 |
IF (rdenom .ne. 0.) THEN |
164 |
Cr = (stencil(0)-stencil(-1))/rdenom |
165 |
ELSE |
166 |
Cr = 1.E20 * (stencil(0)-stencil(-1)) |
167 |
ENDIF |
168 |
|
169 |
IF (STREAMICE_ADV_SCHEME.ne.'DST3') THEN |
170 |
! phi = SLOPE_LIMITER(stencil(0)-stencil(-1), |
171 |
! & stencil(1)-stencil(0)) |
172 |
phi = SLOPE_LIMITER (Cr) |
173 |
|
174 |
ELSE |
175 |
d0 = (2.-cfl)*(1.-cfl)/6.0 |
176 |
d1 = (1.-cfl**2)/6.0 |
177 |
psi = d0+d1*Cr |
178 |
phi = MAX(0. _d 0,MIN(MIN(1. _d 0,psi), |
179 |
& Cr*(1. _d 0 -CFL)/(CFL+1. _d -20) )) |
180 |
ENDIF |
181 |
|
182 |
IF (STREAMICE_ADV_SCHEME.ne.'DST3') THEN |
183 |
XFLUX (i,j,bi,bj) = uface * |
184 |
& (stencil(0) + phi * .5 * (1.0-cfl) * |
185 |
& (stencil(1)-stencil(0))) |
186 |
ELSE |
187 |
XFLUX (i,j,bi,bj) = uface * |
188 |
& (stencil(0) + phi * |
189 |
& (stencil(1)-stencil(0))) |
190 |
ENDIF |
191 |
|
192 |
ELSE ! one of the two cells needed for a HO scheme is missing, use FO scheme |
193 |
|
194 |
Xflux (i,j,bi,bj) = uface * stencil(0) |
195 |
|
196 |
ENDIF |
197 |
|
198 |
ELSE |
199 |
|
200 |
Xflux (i,j,bi,bj) = 0. _d 0 |
201 |
|
202 |
ENDIF |
203 |
|
204 |
ENDIF |
205 |
|
206 |
ENDIF |
207 |
ENDDO |
208 |
ENDIF |
209 |
ENDDO |
210 |
ENDDO |
211 |
ENDDO |
212 |
|
213 |
!C X-FLUXES AT CELL BOUNDARIES CALCULATED; NOW TAKE FLUX DIVERGENCE TO INCREMENT THICKNESS |
214 |
! |
215 |
! DO bj=myByLo(myThid),myByHi(myThid) |
216 |
! DO bi=myBxLo(myThid),myBxHi(myThid) |
217 |
! DO j=1-3,sNy+3 |
218 |
! Gj = (myYGlobalLo-1)+(bj-1)*sNy+j |
219 |
! IF ((Gj .ge. 1) .and. (Gj .le. Ny)) THEN |
220 |
! DO i=1-2,sNx+2 |
221 |
! IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN |
222 |
! h(i,j,bi,bj) = h(i,j,bi,bj) - time_step * |
223 |
! & (hflux_x(i+1,j,bi,bj)*dyG(i+1,j,bi,bj) - |
224 |
! & hflux_x(i,j,bi,bj)*dyG(i,j,bi,bj)) * |
225 |
! & recip_rA (i,j,bi,bj) |
226 |
! ENDIF |
227 |
! ENDDO |
228 |
! ENDIF |
229 |
! ENDDO |
230 |
! ENDDO |
231 |
! ENDDO |
232 |
|
233 |
#endif |
234 |
RETURN |
235 |
END SUBROUTINE STREAMICE_ADV_FLUX_FL_X |