/[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.10 - (show annotations) (download)
Wed Jan 9 21:56:18 2013 UTC (12 years, 6 months ago) by dgoldberg
Branch: MAIN
Changes since 1.9: +70 -5 lines
changes to accept real datasets

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

  ViewVC Help
Powered by ViewVC 1.1.22