/[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.7 - (show annotations) (download)
Thu Oct 4 15:40:16 2012 UTC (12 years, 9 months ago) by dgoldberg
Branch: MAIN
Changes since 1.6: +2 -1 lines
new field and cost function for thickness drift

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

  ViewVC Help
Powered by ViewVC 1.1.22