/[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.9 - (show annotations) (download)
Wed Nov 21 01:04:21 2012 UTC (12 years, 8 months ago) by dgoldberg
Branch: MAIN
Changes since 1.8: +29 -6 lines
call to solve velocity (and initialize ice viscosity) if AD is turned on; velocities overwritten in first timestep

1 C $Header: /u/gcmpack/MITgcm_contrib/dgoldberg/streamice/streamice_init_varia.F,v 1.8 2012/10/10 15:03:10 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 ! C_basal_friction(i,j,bi,bj) = C_basal_fric_const
64 A_glen(i,j,bi,bj) = A_glen_isothermal
65 H_streamice_prev(i,j,bi,bj) = 0. _d 0
66 #ifdef ALLOW_AUTODIFF_TAMC
67 ru_old_si(i,j,bi,bj) = 0. _d 0
68 rv_old_si(i,j,bi,bj) = 0. _d 0
69 zu_old_si(i,j,bi,bj) = 0. _d 0
70 zv_old_si(i,j,bi,bj) = 0. _d 0
71 h_after_uflux_SI(i,j,bi,bj) = 0. _d 0
72 #endif
73 #ifdef USE_ALT_RLOW
74 R_low_si(i,j,bi,bj) = 0. _d 0
75 #endif
76
77 #ifdef STREAMICE_HYBRID_STRESS
78 do k=1,Nr
79 visc_streamice_full(i,j,k,bi,bj) =
80 & eps_glen_min**((1-n_glen)/n_glen)
81 enddo
82 streamice_taubx (i,j,bi,bj) = 0. _d 0
83 streamice_tauby (i,j,bi,bj) = 0. _d 0
84 #endif
85 ENDDO
86 ENDDO
87
88 #ifdef ALLOW_COST_TEST
89 cost_func1_streamice (bi,bj) = 0.0
90 #endif
91
92 ENDDO
93 ENDDO
94
95 DO j = 1-oly, sNy+oly
96 DO i = 1-olx, sNx+olx
97 DO bj = myByLo(myThid), myByHi(myThid)
98 DO bi = myBxLo(myThid), myBxHi(myThid)
99 cc DO k=1,4
100 DO col_x=-1,1
101 DO col_y=-1,1
102 streamice_cg_A1(i,j,bi,bj,col_x,col_y)=0.0
103 streamice_cg_A2(i,j,bi,bj,col_x,col_y)=0.0
104 streamice_cg_A3(i,j,bi,bj,col_x,col_y)=0.0
105 streamice_cg_A4(i,j,bi,bj,col_x,col_y)=0.0
106 ENDDO
107 ENDDO
108 cc ENDDO
109 ENDDO
110 ENDDO
111 ENDDO
112 ENDDO
113
114 C INIT. INTEGER ARRAYS
115
116 DO bj = myByLo(myThid), myByHi(myThid)
117 DO bi = myBxLo(myThid), myBxHi(myThid)
118 DO j=1-Oly,sNy+Oly
119 DO i=1-Olx,sNx+Olx
120 STREAMICE_hmask(i,j,bi,bj) = -1.0
121 STREAMICE_umask(i,j,bi,bj) = 0.0
122 STREAMICE_vmask(i,j,bi,bj) = 0.0
123 STREAMICE_ufacemask(i,j,bi,bj) = 0.0
124 STREAMICE_vfacemask(i,j,bi,bj) = 0.0
125 STREAMICE_float_cond(i,j,bi,bj) = 0.0
126 ENDDO
127 ENDDO
128 ENDDO
129 ENDDO
130
131
132 #ifdef USE_ALT_RLOW
133 ! init alternate array for topog
134 IF ( bathyFile .NE. ' ' ) THEN
135 _BARRIER
136 C The 0 is the "iteration" argument. The ' ' is an empty suffix
137 CALL READ_FLD_XY_RS( bathyFile, '',
138 & R_low_si, 0, myThid )
139
140 ELSE
141 WRITE(msgBuf,'(A)') 'STREAMICE TOPOG - FILENAME MISSING'
142 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
143 & SQUEEZE_RIGHT , 1)
144 ENDIF
145 #endif
146
147 ! initialize thickness
148
149 IF ( STREAMICEthickInit.EQ.'PARAM' ) THEN
150
151 WRITE(msgBuf,'(A)') 'initializing analytic thickness'
152 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
153 & SQUEEZE_RIGHT , 1)
154
155 slope_pos = shelf_edge_pos - shelf_flat_width
156 c1 = 0.0
157 IF (shelf_slope_scale .GT. 0.0) THEN
158 c1 = 1.0 / shelf_slope_scale
159 ENDIF
160
161 DO bj = myByLo(myThid), myByHi(myThid)
162 DO bi = myBxLo(myThid), myBxHi(myThid)
163 DO j=1,sNy
164 DO i=1,sNx
165 Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
166 Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
167
168 IF ((Gi.lt.Nx).and.(Gj.lt.Ny)) THEN
169
170 C IF (flow_dir .EQ. 2.0) THEN
171 IF (.TRUE.) THEN
172 IF (xC(i-1,j,bi,bj).GE.shelf_edge_pos) THEN
173 area_shelf_streamice(i,j,bi,bj) = 0. _d 0
174 STREAMICE_hmask(i,j,bi,bj) = 0. _d 0
175 ELSE
176
177 IF (xC(i,j,bi,bj).GT.slope_pos) THEN
178 H_streamice (i,j,bi,bj) = shelf_min_draft
179 ELSE
180 H_streamice (i,j,bi,bj) = (shelf_min_draft +
181 & (shelf_max_draft - shelf_min_draft) *
182 & min (1.0, (c1*(slope_pos-xC(i,j,bi,bj)))**2))
183 ENDIF
184
185 IF (xC(i,j,bi,bj).GT.shelf_edge_pos) THEN
186 area_shelf_streamice(i,j,bi,bj) = rA(i,j,bi,bj) *
187 & (shelf_edge_pos-xG(i,j,bi,bj)) /
188 & (xG(i+1,j,bi,bj)-xG(i,j,bi,bj))
189 IF (area_shelf_streamice(i,j,bi,bj).gt. 0. _d 0) THEN
190 STREAMICE_hmask(i,j,bi,bj) = 2.0
191 ELSE
192 STREAMICE_hmask(i,j,bi,bj) = 0.0
193 H_streamice(i,j,bi,bj) = 0.0
194 ENDIF
195 ELSE
196 area_shelf_streamice(i,j,bi,bj) = rA(i,j,bi,bj)
197 STREAMICE_hmask(i,j,bi,bj) = 1.0
198 ENDIF
199
200
201 ENDIF
202 ENDIF
203 ENDIF
204 ENDDO
205 ENDDO
206 ENDDO
207 ENDDO
208
209 ELSE IF ( STREAMICEthickInit.EQ.'FILE' ) THEN
210
211 IF ( STREAMICEthickFile .NE. ' ' ) THEN
212 _BARRIER
213 C The 0 is the "iteration" argument. The ' ' is an empty suffix
214 CALL READ_FLD_XY_RS( STREAMICEthickFile, ' ', H_streamice,
215 & 0, myThid )
216 DO bj = myByLo(myThid), myByHi(myThid)
217 DO bi = myBxLo(myThid), myBxHi(myThid)
218 DO j=1,sNy
219 DO i=1,sNx
220 Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
221 Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
222 IF ((Gi.lt.Nx.OR.STREAMICE_EW_periodic).and.
223 & (Gj.lt.Ny.OR.STREAMICE_NS_periodic)) THEN
224 IF (H_streamice(i,j,bi,bj).GT.0. _d 0) THEN
225 area_shelf_streamice(i,j,bi,bj) = rA(i,j,bi,bj)
226 STREAMICE_hmask(i,j,bi,bj) = 1.0
227 ELSE
228 area_shelf_streamice(i,j,bi,bj) = 0. _d 0
229 STREAMICE_hmask(i,j,bi,bj) = 0. _d 0
230 ENDIF
231 Do k=1,Nr
232 STREAMICE_ctrl_mask(i,j,bi,bj,k) = 1. _d 0
233 enddo
234 ENDIF
235 ENDDO
236 ENDDO
237 ENDDO
238 ENDDO
239 ELSE
240 WRITE(msgBuf,'(A)') 'INIT THICKNESS - FILENAME MISSING'
241 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
242 & SQUEEZE_RIGHT , 1)
243 ENDIF
244
245 ELSE
246
247 WRITE(msgBuf,'(A)') 'INIT THICKNESS - NOT IMPLENTED'
248 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
249 & SQUEEZE_RIGHT , 1)
250 ENDIF
251
252 ! finish initialize thickness
253
254 ! initialize basal traction
255
256 IF ( STREAMICEbasalTracConfig.EQ.'FILE' ) THEN
257
258 IF ( STREAMICEbasalTracFile .NE. ' ' ) THEN
259 _BARRIER
260 C The 0 is the "iteration" argument. The ' ' is an empty suffix
261 CALL READ_FLD_XY_RS( STREAMICEbasalTracFile, ' ',
262 & C_basal_friction, 0, myThid )
263
264 ELSE
265 WRITE(msgBuf,'(A)') 'INIT THICKNESS - FILENAME MISSING'
266 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
267 & SQUEEZE_RIGHT , 1)
268 ENDIF
269
270 ELSE IF (STREAMICEbasalTracConfig.EQ.'UNIFORM' ) THEN
271
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 C_basal_friction(i,j,bi,bj) = C_basal_fric_const
277 ENDDO
278 ENDDO
279 ENDDO
280 ENDDO
281
282 ELSE IF (STREAMICEbasalTracConfig.EQ.'2DPERIODIC' ) THEN
283
284 lenx = sNx*nSx*nPx*delX(1)
285 leny = sNy*nSy*nPy*delY(1)
286 print *, 'lenx', lenx
287 print *, 'leny', leny
288 DO bj = myByLo(myThid), myByHi(myThid)
289 DO bi = myBxLo(myThid), myBxHi(myThid)
290 DO j=1,sNy
291 DO i=1,sNx
292 x = xC(i,j,bi,bj)
293 y = yC(i,j,bi,bj)
294 C_basal_friction(i,j,bi,bj) =
295 & sqrt(C_basal_fric_const**2*
296 & (1+sin(2*streamice_kx_b_init*PI*x/lenx)*
297 & sin(2*streamice_ky_b_init*PI*y/leny)))
298 ENDDO
299 ENDDO
300 ENDDO
301 ENDDO
302
303 ELSE IF (STREAMICEbasalTracConfig.EQ.'1DPERIODIC' ) THEN
304
305 lenx = sNx*nSx*nPx*delX(1)
306 print *, 'lenx', lenx
307 DO bj = myByLo(myThid), myByHi(myThid)
308 DO bi = myBxLo(myThid), myBxHi(myThid)
309 DO j=1,sNy
310 DO i=1,sNx
311 x = xC(i,j,bi,bj)
312 y = yC(i,j,bi,bj)
313 C_basal_friction(i,j,bi,bj) =
314 & sqrt(C_basal_fric_const**2*(1+
315 & sin(2*streamice_kx_b_init*PI*x/lenx)))
316 ENDDO
317 ENDDO
318 ENDDO
319 ENDDO
320
321 ELSE
322
323 WRITE(msgBuf,'(A)') 'INIT TRAC - NOT IMPLENTED'
324 CALL PRINT_ERROR( msgBuf, myThid)
325 STOP 'ABNORMAL END: S/R STREAMICE_INIT_VAR'
326 ENDIF
327
328 ! finish initialize basal traction
329
330 CALL STREAMICE_UPD_FFRAC_UNCOUPLED ( myThid )
331
332 _EXCH_XY_RL(H_streamice, myThid )
333 _EXCH_XY_RL(STREAMICE_hmask, myThid )
334 _EXCH_XY_RL(area_shelf_streamice, myThid )
335 _EXCH_XY_RL(C_basal_friction, myThid )
336 #ifdef USE_ALT_RLOW
337 _EXCH_XY_RL(R_low_si, myThid )
338 #endif
339
340 !#ifdef STREAMICE_HYBRID_STRESS
341
342 ! CALL STREAMICE_VISC_BETA (myThid)
343
344 ! DNG THIS CALL IS TO INITIALISE VISCOSITY
345 ! TO AVOID POSSIBLE ADJOINT INSTABILITIES
346 ! IT IS WRITTEN OVER IN FIRST TIMESTEP
347
348 #ifdef ALLOW_AUTODIFF
349
350 CALL STREAMICE_UPD_FFRAC_UNCOUPLED ( myThid )
351 CALL STREAMICE_VELMASK_UPD (myThid)
352 CALL STREAMICE_VEL_SOLVE( myThid )
353
354 #endif
355
356 !#endif
357
358 CALL WRITE_FLD_XY_RL ( "C_basal_fric", "",
359 & C_basal_friction, 0, myThid )
360 CALL WRITE_FLD_XY_RL ( "H_streamice", "init",
361 & H_streamIce, 0, myThid )
362 CALL WRITE_FLD_XY_RL ( "area_shelf_streamice", "init",
363 & area_shelf_streamice, 0, myThid )
364 CALL WRITE_FLD_XY_RL ( "STREAMICE_hmask", "init",
365 & STREAMICE_hmask, 0, myThid )
366 #ifdef ALLOW_CTRL
367 CALL ACTIVE_WRITE_GEN_RS( 'maskCtrlst', STREAMICE_ctrl_mask,
368 & 'XY', Nr, 1, .FALSE., 0, mythid, dummyRS )
369 #endif
370 ! call active_write_xyz( 'maskCtrlS', STREAMICE_ctrl_mask, 1, 0,
371 ! & mythid, dummy)
372 ! CALL STREAMICE_VELMASK_UPD (myThid)
373 ! CALL STREAMICE_UPD_FFRAC_UNCOUPLED ( myThid )
374 ! CALL STREAMICE_VEL_SOLVE( myThid )
375
376 CALL WRITE_FLD_XY_RL ( "U_init", "",
377 & C_basal_friction, 0, myThid )
378 CALL WRITE_FLD_XY_RL ( "V_init", "",
379 & V_streamice, 0, myThid )
380 #ifdef USE_ALT_RLOW
381 CALL WRITE_FLD_XY_RL ( "R_low_si", "init",
382 & R_low_si, 0, myThid )
383 #endif
384
385 ! CALL WRITE_FULLARRAY_RL ("H",H_streamice,1,0,0,1,0,myThid)
386 ! CALL WRITE_FULLARRAY_RL ("hmask",STREAMICE_hmask,1,0,0,1,0,myThid)
387 ! CALL WRITE_FULLARRAY_RL ("umask",STREAMICE_umask,1,0,0,1,0,myThid)
388
389 #endif /* ALLOW_STREAMICE */
390
391 RETURN
392 END
393

  ViewVC Help
Powered by ViewVC 1.1.22