/[MITgcm]/MITgcm_contrib/dgoldberg/streamice/streamice_init_varia.F
ViewVC logotype

Contents of /MITgcm_contrib/dgoldberg/streamice/streamice_init_varia.F

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


Revision 1.4 - (show annotations) (download)
Tue Sep 18 17:06:48 2012 UTC (12 years, 10 months ago) by dgoldberg
Branch: MAIN
Changes since 1.3: +101 -9 lines
changes for periodic boundary conds and hybrid stress balance

1 C $Header: /u/gcmpack/MITgcm_contrib/dgoldberg/streamice/streamice_init_varia.F,v 1.3 2012/09/04 21:11:44 dgoldberg Exp $
2 C $Name: $
3
4 #include "STREAMICE_OPTIONS.h"
5
6 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
7
8 CBOP
9 SUBROUTINE STREAMICE_INIT_VARIA( myThid )
10 C /============================================================\
11 C | SUBROUTINE STREAMICE_INIT_VARIA |
12 C | o Routine to initialize STREAMICE variables. |
13 C |============================================================|
14 C | Initialize STREAMICE parameters and variables. |
15 C \============================================================/
16 IMPLICIT NONE
17
18 C === Global variables ===
19 #include "SIZE.h"
20 #include "GRID.h"
21 #include "SET_GRID.h"
22 #include "EEPARAMS.h"
23 #include "PARAMS.h"
24 #include "STREAMICE.h"
25 #include "STREAMICE_CG.h"
26 #include "STREAMICE_ADV.h"
27
28 C === Routine arguments ===
29 C myThid - Number of this instance of STREAMICE_INIT_VARIA
30 INTEGER myThid
31 CEndOfInterface
32
33 #ifdef ALLOW_STREAMICE
34 C === Local variables ===
35 C I,J,bi,bj - Loop counters
36 INTEGER i, j, k, bi, bj, Gi, Gj
37 INTEGER col_y, col_x
38 _RL slope_pos, c1, x, y, lenx, leny
39 CHARACTER*(MAX_LEN_MBUF) msgBuf
40 CEOP
41
42 C ZERO OUT FLOATING POINT ARRAYS
43
44 DO bj = myByLo(myThid), myByHi(myThid)
45 DO bi = myBxLo(myThid), myBxHi(myThid)
46 DO j=1-Oly,sNy+Oly
47 DO i=1-Olx,sNx+Olx
48 H_streamIce(i,j,bi,bj) = 0. _d 0
49 U_streamice(i,j,bi,bj) = 0. _d 0
50 V_streamice(i,j,bi,bj) = 0. _d 0
51 visc_streamice(i,j,bi,bj) = 0. _d 0
52 tau_beta_eff_streamice(i,j,bi,bj) = 0. _d 0
53 float_frac_streamice(i,j,bi,bj) = 0. _d 0
54 base_el_streamice(i,j,bi,bj) = 0. _d 0
55 surf_el_streamice(i,j,bi,bj) = 0. _d 0
56 area_shelf_streamice(i,j,bi,bj) = 0. _d 0
57 mass_ice_streamice(i,j,bi,bj) = 0. _d 0
58 BDOT_streamice(i,j,bi,bj) = 0. _d 0
59 ! C_basal_friction(i,j,bi,bj) = C_basal_fric_const
60 A_glen(i,j,bi,bj) = A_glen_isothermal
61 #ifdef ALLOW_AUTODIFF_TAMC
62 ru_old_si(i,j,bi,bj) = 0. _d 0
63 rv_old_si(i,j,bi,bj) = 0. _d 0
64 zu_old_si(i,j,bi,bj) = 0. _d 0
65 zv_old_si(i,j,bi,bj) = 0. _d 0
66 h_after_uflux_SI(i,j,bi,bj) = 0. _d 0
67 #endif
68 #ifdef STREAMICE_HYBRID_STRESS
69 do k=1,Nr
70 visc_streamice_full(i,j,k,bi,bj) =
71 & eps_glen_min**((1-n_glen)/n_glen)
72 enddo
73 #endif
74 ENDDO
75 ENDDO
76 ENDDO
77 ENDDO
78
79 DO j = 1-oly, sNy+oly
80 DO i = 1-olx, sNx+olx
81 DO bj = myByLo(myThid), myByHi(myThid)
82 DO bi = myBxLo(myThid), myBxHi(myThid)
83 cc DO k=1,4
84 DO col_x=-1,1
85 DO col_y=-1,1
86 streamice_cg_A1(i,j,bi,bj,col_x,col_y)=0.0
87 streamice_cg_A2(i,j,bi,bj,col_x,col_y)=0.0
88 streamice_cg_A3(i,j,bi,bj,col_x,col_y)=0.0
89 streamice_cg_A4(i,j,bi,bj,col_x,col_y)=0.0
90 ENDDO
91 ENDDO
92 cc ENDDO
93 ENDDO
94 ENDDO
95 ENDDO
96 ENDDO
97
98 C INIT. INTEGER ARRAYS
99
100 DO bj = myByLo(myThid), myByHi(myThid)
101 DO bi = myBxLo(myThid), myBxHi(myThid)
102 DO j=1-Oly,sNy+Oly
103 DO i=1-Olx,sNx+Olx
104 STREAMICE_hmask(i,j,bi,bj) = -1.0
105 STREAMICE_umask(i,j,bi,bj) = 0.0
106 STREAMICE_vmask(i,j,bi,bj) = 0.0
107 STREAMICE_ufacemask(i,j,bi,bj) = 0.0
108 STREAMICE_vfacemask(i,j,bi,bj) = 0.0
109 STREAMICE_float_cond(i,j,bi,bj) = 0.0
110 ENDDO
111 ENDDO
112 ENDDO
113 ENDDO
114
115 ! initialize thickness
116
117 IF ( STREAMICEthickInit.EQ.'PARAM' ) THEN
118
119 WRITE(msgBuf,'(A)') 'initializing analytic thickness'
120 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
121 & SQUEEZE_RIGHT , 1)
122
123 slope_pos = shelf_edge_pos - shelf_flat_width
124 c1 = 0.0
125 IF (shelf_slope_scale .GT. 0.0) THEN
126 c1 = 1.0 / shelf_slope_scale
127 ENDIF
128
129 DO bj = myByLo(myThid), myByHi(myThid)
130 DO bi = myBxLo(myThid), myBxHi(myThid)
131 DO j=1,sNy
132 DO i=1,sNx
133 Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
134 Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
135
136 IF ((Gi.lt.Nx).and.(Gj.lt.Ny)) THEN
137
138 C IF (flow_dir .EQ. 2.0) THEN
139 IF (.TRUE.) THEN
140 IF (xC(i-1,j,bi,bj).GE.shelf_edge_pos) THEN
141 area_shelf_streamice(i,j,bi,bj) = 0. _d 0
142 STREAMICE_hmask(i,j,bi,bj) = 0. _d 0
143 ELSE
144
145 IF (xC(i,j,bi,bj).GT.slope_pos) THEN
146 H_streamice (i,j,bi,bj) = shelf_min_draft
147 ELSE
148 H_streamice (i,j,bi,bj) = (shelf_min_draft +
149 & (shelf_max_draft - shelf_min_draft) *
150 & min (1.0, (c1*(slope_pos-xC(i,j,bi,bj)))**2))
151 ENDIF
152
153 IF (xC(i,j,bi,bj).GT.shelf_edge_pos) THEN
154 area_shelf_streamice(i,j,bi,bj) = rA(i,j,bi,bj) *
155 & (shelf_edge_pos-xG(i,j,bi,bj)) /
156 & (xG(i+1,j,bi,bj)-xG(i,j,bi,bj))
157 IF (area_shelf_streamice(i,j,bi,bj).gt. 0._d 0) THEN
158 STREAMICE_hmask(i,j,bi,bj) = 2.0
159 ELSE
160 STREAMICE_hmask(i,j,bi,bj) = 0.0
161 H_streamice(i,j,bi,bj) = 0.0
162 ENDIF
163 ELSE
164 area_shelf_streamice(i,j,bi,bj) = rA(i,j,bi,bj)
165 STREAMICE_hmask(i,j,bi,bj) = 1.0
166 ENDIF
167
168
169 ENDIF
170 ENDIF
171 ENDIF
172 ENDDO
173 ENDDO
174 ENDDO
175 ENDDO
176
177 ELSE IF ( STREAMICEthickInit.EQ.'FILE' ) THEN
178
179 IF ( STREAMICEthickFile .NE. ' ' ) THEN
180 _BARRIER
181 C The 0 is the "iteration" argument. The ' ' is an empty suffix
182 CALL READ_FLD_XY_RS( STREAMICEthickFile, ' ', H_streamice,
183 & 0, myThid )
184 DO bj = myByLo(myThid), myByHi(myThid)
185 DO bi = myBxLo(myThid), myBxHi(myThid)
186 DO j=1,sNy
187 DO i=1,sNx
188 Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
189 Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
190 IF ((Gi.lt.Nx.OR.STREAMICE_EW_periodic).and.
191 & (Gj.lt.Ny.OR.STREAMICE_NS_periodic)) THEN
192 IF (H_streamice(i,j,bi,bj).GT.0. _d 0) THEN
193 area_shelf_streamice(i,j,bi,bj) = rA(i,j,bi,bj)
194 STREAMICE_hmask(i,j,bi,bj) = 1.0
195 ELSE
196 area_shelf_streamice(i,j,bi,bj) = 0. _d 0
197 STREAMICE_hmask(i,j,bi,bj) = 0. _d 0
198 ENDIF
199 ENDIF
200 ENDDO
201 ENDDO
202 ENDDO
203 ENDDO
204 ELSE
205 WRITE(msgBuf,'(A)') 'INIT THICKNESS - FILENAME MISSING'
206 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
207 & SQUEEZE_RIGHT , 1)
208 ENDIF
209
210 ELSE
211
212 WRITE(msgBuf,'(A)') 'INIT THICKNESS - NOT IMPLENTED'
213 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
214 & SQUEEZE_RIGHT , 1)
215 ENDIF
216
217 ! finish initialize thickness
218
219 ! initialize basal traction
220
221 IF ( STREAMICEbasalTracConfig.EQ.'FILE' ) THEN
222
223 IF ( STREAMICEbasalTracFile .NE. ' ' ) THEN
224 _BARRIER
225 C The 0 is the "iteration" argument. The ' ' is an empty suffix
226 CALL READ_FLD_XY_RS( STREAMICEbasalTracFile, ' ',
227 & C_basal_friction, 0, myThid )
228
229 ELSE
230 WRITE(msgBuf,'(A)') 'INIT THICKNESS - FILENAME MISSING'
231 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
232 & SQUEEZE_RIGHT , 1)
233 ENDIF
234
235 ELSE IF (STREAMICEbasalTracConfig.EQ.'UNIFORM' ) THEN
236
237 DO bj = myByLo(myThid), myByHi(myThid)
238 DO bi = myBxLo(myThid), myBxHi(myThid)
239 DO j=1,sNy
240 DO i=1,sNx
241 C_basal_friction(i,j,bi,bj) = C_basal_fric_const
242 ENDDO
243 ENDDO
244 ENDDO
245 ENDDO
246
247 ELSE IF (STREAMICEbasalTracConfig.EQ.'2DPERIODIC' ) THEN
248
249 lenx = sNx*nSx*nPx*delX(1)
250 leny = sNy*nSy*nPy*delY(1)
251 print *, 'lenx', lenx
252 print *, 'leny', leny
253 DO bj = myByLo(myThid), myByHi(myThid)
254 DO bi = myBxLo(myThid), myBxHi(myThid)
255 DO j=1,sNy
256 DO i=1,sNx
257 x = xC(i,j,bi,bj)
258 y = yC(i,j,bi,bj)
259 C_basal_friction(i,j,bi,bj) =
260 & sqrt(C_basal_fric_const**2*
261 & (1+sin(2*streamice_kx_b_init*PI*x/lenx)*
262 & sin(2*streamice_ky_b_init*PI*y/leny)))
263 ENDDO
264 ENDDO
265 ENDDO
266 ENDDO
267
268 ELSE IF (STREAMICEbasalTracConfig.EQ.'1DPERIODIC' ) THEN
269
270 lenx = sNx*nSx*nPx*delX(1)
271 print *, 'lenx', lenx
272 DO bj = myByLo(myThid), myByHi(myThid)
273 DO bi = myBxLo(myThid), myBxHi(myThid)
274 DO j=1,sNy
275 DO i=1,sNx
276 x = xC(i,j,bi,bj)
277 y = yC(i,j,bi,bj)
278 C_basal_friction(i,j,bi,bj) =
279 & sqrt(C_basal_fric_const**2*(1+
280 & sin(2*streamice_kx_b_init*PI*x/lenx)))
281 ENDDO
282 ENDDO
283 ENDDO
284 ENDDO
285
286 ELSE
287
288 WRITE(msgBuf,'(A)') 'INIT THICKNESS - NOT IMPLENTED'
289 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
290 & SQUEEZE_RIGHT , 1)
291 ENDIF
292
293 ! finish initialize basal traction
294
295 CALL STREAMICE_UPD_FFRAC_UNCOUPLED ( myThid )
296
297 _EXCH_XY_RL(H_streamice, myThid )
298 _EXCH_XY_RL(STREAMICE_hmask, myThid )
299 _EXCH_XY_RL(area_shelf_streamice, myThid )
300 _EXCH_XY_RL(C_basal_friction, myThid )
301
302
303 #ifdef STREAMICE_HYBRID_STRESS
304
305 CALL STREAMICE_VISC_BETA (myThid)
306
307 #endif
308
309 CALL WRITE_FLD_XY_RL ( "H_streamice", "init",
310 & H_streamIce, 0, myThid )
311 CALL WRITE_FLD_XY_RL ( "area_shelf_streamice", "init",
312 & area_shelf_streamice, 0, myThid )
313 CALL WRITE_FLD_XY_RL ( "STREAMICE_hmask", "init",
314 & STREAMICE_hmask, 0, myThid )
315
316 ! CALL STREAMICE_VELMASK_UPD (myThid)
317 ! CALL STREAMICE_UPD_FFRAC_UNCOUPLED ( myThid )
318 ! CALL STREAMICE_VEL_SOLVE( myThid )
319
320 CALL WRITE_FLD_XY_RL ( "U_init", "",
321 & C_basal_friction, 0, myThid )
322 CALL WRITE_FLD_XY_RL ( "V_init", "",
323 & V_streamice, 0, myThid )
324
325 ! CALL WRITE_FULLARRAY_RL ("H",H_streamice,1,0,0,1,0,myThid)
326 ! CALL WRITE_FULLARRAY_RL ("hmask",STREAMICE_hmask,1,0,0,1,0,myThid)
327 ! CALL WRITE_FULLARRAY_RL ("umask",STREAMICE_umask,1,0,0,1,0,myThid)
328
329 #endif /* ALLOW_STREAMICE */
330
331 RETURN
332 END
333

  ViewVC Help
Powered by ViewVC 1.1.22