/[MITgcm]/MITgcm/pkg/streamice/streamice_init_varia.F
ViewVC logotype

Contents of /MITgcm/pkg/streamice/streamice_init_varia.F

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


Revision 1.3 - (show annotations) (download)
Fri Jun 21 21:15:36 2013 UTC (10 years, 11 months ago) by dgoldberg
Branch: MAIN
Changes since 1.2: +2 -2 lines
replace EXCH_XY_RL by EXCH_XY_RS for RS arrays

1 C $Header: /u/gcmpack/MITgcm/pkg/streamice/streamice_init_varia.F,v 1.2 2013/06/21 17:53:25 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
49 DO bj = myByLo(myThid), myByHi(myThid)
50 DO bi = myBxLo(myThid), myBxHi(myThid)
51 DO j=1-Oly,sNy+Oly
52 DO i=1-Olx,sNx+Olx
53 H_streamIce(i,j,bi,bj) = 0. _d 0
54 U_streamice(i,j,bi,bj) = 0. _d 0
55 V_streamice(i,j,bi,bj) = 0. _d 0
56 visc_streamice(i,j,bi,bj) = 0. _d 0
57 tau_beta_eff_streamice(i,j,bi,bj) = 0. _d 0
58 float_frac_streamice(i,j,bi,bj) = 0. _d 0
59 base_el_streamice(i,j,bi,bj) = 0. _d 0
60 surf_el_streamice(i,j,bi,bj) = 0. _d 0
61 area_shelf_streamice(i,j,bi,bj) = 0. _d 0
62 mass_ice_streamice(i,j,bi,bj) = 0. _d 0
63 BDOT_streamice(i,j,bi,bj) = 0. _d 0
64 ADOT_streamice(i,j,bi,bj) = streamice_adot_uniform
65 C_basal_friction(i,j,bi,bj) = C_basal_fric_const
66 B_glen(i,j,bi,bj) = B_glen_isothermal
67 H_streamice_prev(i,j,bi,bj) = 0. _d 0
68 #ifdef ALLOW_STREAMICE_2DTRACER
69 #ifdef STREAMICE_TRACER_AB
70 GAD_trac_2d (i,j,bi,bj) = 0. _d 0
71 ! GAD_trac_2dNm1 (i,j,bi,bj) = 0. _d 0
72 #endif
73 #endif
74 #ifdef ALLOW_AUTODIFF_TAMC
75 ru_old_si(i,j,bi,bj) = 0. _d 0
76 rv_old_si(i,j,bi,bj) = 0. _d 0
77 zu_old_si(i,j,bi,bj) = 0. _d 0
78 zv_old_si(i,j,bi,bj) = 0. _d 0
79 ! h_after_uflux_SI(i,j,bi,bj) = 0. _d 0
80 #endif
81 #ifdef USE_ALT_RLOW
82 R_low_si(i,j,bi,bj) = 0. _d 0
83 #endif
84
85 #ifdef STREAMICE_HYBRID_STRESS
86 do k=1,Nr
87 visc_streamice_full(i,j,k,bi,bj) =
88 & eps_glen_min**((1-n_glen)/n_glen)
89 enddo
90 streamice_taubx (i,j,bi,bj) = 0. _d 0
91 streamice_tauby (i,j,bi,bj) = 0. _d 0
92 #endif
93 ENDDO
94 ENDDO
95
96 #ifdef ALLOW_COST_TEST
97 cost_func1_streamice (bi,bj) = 0.0
98 #endif
99
100 ENDDO
101 ENDDO
102
103 DO j = 1-oly, sNy+oly
104 DO i = 1-olx, sNx+olx
105 DO bj = myByLo(myThid), myByHi(myThid)
106 DO bi = myBxLo(myThid), myBxHi(myThid)
107 cc DO k=1,4
108 DO col_x=-1,1
109 DO col_y=-1,1
110 streamice_cg_A1(i,j,bi,bj,col_x,col_y)=0.0
111 streamice_cg_A2(i,j,bi,bj,col_x,col_y)=0.0
112 streamice_cg_A3(i,j,bi,bj,col_x,col_y)=0.0
113 streamice_cg_A4(i,j,bi,bj,col_x,col_y)=0.0
114 ENDDO
115 ENDDO
116 cc ENDDO
117 ENDDO
118 ENDDO
119 ENDDO
120 ENDDO
121
122 C INIT. INTEGER ARRAYS
123
124 DO bj = myByLo(myThid), myByHi(myThid)
125 DO bi = myBxLo(myThid), myBxHi(myThid)
126 DO j=1-Oly,sNy+Oly
127 DO i=1-Olx,sNx+Olx
128 STREAMICE_hmask(i,j,bi,bj) = -1.0
129 STREAMICE_umask(i,j,bi,bj) = 0.0
130 STREAMICE_vmask(i,j,bi,bj) = 0.0
131 STREAMICE_ufacemask(i,j,bi,bj) = 0.0
132 STREAMICE_vfacemask(i,j,bi,bj) = 0.0
133 STREAMICE_float_cond(i,j,bi,bj) = 0.0
134 ENDDO
135 ENDDO
136 ENDDO
137 ENDDO
138
139
140 #ifdef USE_ALT_RLOW
141 ! init alternate array for topog
142 IF ( STREAMICEtopogFile .NE. ' ' ) THEN
143 _BARRIER
144 C The 0 is the "iteration" argument. The ' ' is an empty suffix
145 CALL READ_FLD_XY_RS( STREAMICEtopogFile, '',
146 & R_low_si, 0, myThid )
147
148 ELSE
149 WRITE(msgBuf,'(A)') 'STREAMICE TOPOG - FILENAME MISSING'
150 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
151 & SQUEEZE_RIGHT , 1)
152 ENDIF
153 #endif
154
155 ! initialize thickness
156
157 #ifndef STREAMICE_GEOM_FILE_SETUP
158
159 IF ( STREAMICEthickInit.EQ.'PARAM' ) THEN
160
161 WRITE(msgBuf,'(A)') 'initializing analytic thickness'
162 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
163 & SQUEEZE_RIGHT , 1)
164
165 slope_pos = shelf_edge_pos - shelf_flat_width
166 c1 = 0.0
167 IF (shelf_slope_scale .GT. 0.0) THEN
168 c1 = 1.0 / shelf_slope_scale
169 ENDIF
170
171 DO bj = myByLo(myThid), myByHi(myThid)
172 DO bi = myBxLo(myThid), myBxHi(myThid)
173 DO j=1,sNy
174 DO i=1,sNx
175 Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
176 Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
177
178 IF ((Gi.lt.Nx).and.(Gj.lt.Ny)) THEN
179
180 C IF (flow_dir .EQ. 2.0) THEN
181 IF (.TRUE.) THEN
182 IF (xC(i-1,j,bi,bj).GE.shelf_edge_pos) THEN
183 area_shelf_streamice(i,j,bi,bj) = 0. _d 0
184 STREAMICE_hmask(i,j,bi,bj) = 0. _d 0
185 ELSE
186
187 IF (xC(i,j,bi,bj).GT.slope_pos) THEN
188 H_streamice (i,j,bi,bj) = shelf_min_draft
189 ELSE
190 H_streamice (i,j,bi,bj) = (shelf_min_draft +
191 & (shelf_max_draft - shelf_min_draft) *
192 & min (oneRL, (c1*(slope_pos-xC(i,j,bi,bj)))**2))
193 ENDIF
194
195 IF (xC(i,j,bi,bj).GT.shelf_edge_pos) THEN
196 area_shelf_streamice(i,j,bi,bj) = rA(i,j,bi,bj) *
197 & (shelf_edge_pos-xG(i,j,bi,bj)) /
198 & (xG(i+1,j,bi,bj)-xG(i,j,bi,bj))
199 IF (area_shelf_streamice(i,j,bi,bj).gt. 0. _d 0) THEN
200 STREAMICE_hmask(i,j,bi,bj) = 2.0
201 ELSE
202 STREAMICE_hmask(i,j,bi,bj) = 0.0
203 H_streamice(i,j,bi,bj) = 0.0
204 ENDIF
205 ELSE
206 area_shelf_streamice(i,j,bi,bj) = rA(i,j,bi,bj)
207 STREAMICE_hmask(i,j,bi,bj) = 1.0
208 ENDIF
209
210
211 ENDIF
212 ENDIF
213 ENDIF
214 ENDDO
215 ENDDO
216 ENDDO
217 ENDDO
218
219 ELSE IF ( STREAMICEthickInit.EQ.'FILE' ) THEN
220
221 IF ( STREAMICEthickFile .NE. ' ' ) THEN
222 _BARRIER
223 C The 0 is the "iteration" argument. The ' ' is an empty suffix
224 CALL READ_FLD_XY_RL( STREAMICEthickFile, ' ', H_streamice,
225 & 0, myThid )
226 DO bj = myByLo(myThid), myByHi(myThid)
227 DO bi = myBxLo(myThid), myBxHi(myThid)
228 DO j=1,sNy
229 DO i=1,sNx
230 Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
231 Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
232 IF ((Gi.lt.Nx.OR.STREAMICE_EW_periodic).and.
233 & (Gj.lt.Ny.OR.STREAMICE_NS_periodic)) THEN
234 IF (H_streamice(i,j,bi,bj).GT.0. _d 0) THEN
235 area_shelf_streamice(i,j,bi,bj) = rA(i,j,bi,bj)
236 STREAMICE_hmask(i,j,bi,bj) = 1.0
237 ELSE
238 area_shelf_streamice(i,j,bi,bj) = 0. _d 0
239 STREAMICE_hmask(i,j,bi,bj) = 0. _d 0
240 ENDIF
241 Do k=1,Nr
242 STREAMICE_ctrl_mask(i,j,k,bi,bj) = 1. _d 0
243 enddo
244 ENDIF
245 ENDDO
246 ENDDO
247 ENDDO
248 ENDDO
249 ELSE
250 WRITE(msgBuf,'(A)') 'INIT THICKNESS - FILENAME MISSING'
251 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
252 & SQUEEZE_RIGHT , 1)
253 ENDIF
254
255 ELSE
256
257 WRITE(msgBuf,'(A)') 'INIT THICKNESS - NOT IMPLENTED'
258 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
259 & SQUEEZE_RIGHT , 1)
260 ENDIF
261
262 #else
263 ! STREAMICE_GEOM_FILE_SETUP - init thickness and hmask MUST come from file
264
265 IF ( STREAMICEthickFile .NE. ' ' ) THEN
266 _BARRIER
267 C The 0 is the "iteration" argument. The ' ' is an empty suffix
268 CALL READ_FLD_XY_RL( STREAMICEthickFile, ' ', H_streamice,
269 & 0, myThid )
270 ELSE
271 WRITE(msgBuf,'(A)') 'INIT THICKNESS - FILENAME MISSING'
272 CALL PRINT_ERROR( msgBuf, myThid)
273 ENDIF
274
275 IF ( STREAMICEhMaskFile .NE. ' ' ) THEN
276 _BARRIER
277 C The 0 is the "iteration" argument. The ' ' is an empty suffix
278 CALL READ_FLD_XY_RS( STREAMICEhMaskFile, ' ', STREAMICE_hmask,
279 & 0, myThid )
280 ELSE
281 WRITE(msgBuf,'(A)') 'INIT HMASK - FILENAME MISSING'
282 CALL PRINT_ERROR( msgBuf, myThid)
283 ENDIF
284
285 #endif
286 ! STREAMICE_GEOM_FILE_SETUP
287
288
289 ! finish initialize thickness
290
291 ! initialize glen constant
292
293 IF ( STREAMICEGlenConstConfig.EQ.'FILE' ) THEN
294
295 IF ( STREAMICEGlenConstFile .NE. ' ' ) THEN
296 _BARRIER
297
298 CALL READ_FLD_XY_RL( STREAMICEGlenConstFile, ' ',
299 & B_glen, 0, myThid )
300
301 ELSE
302 WRITE(msgBuf,'(A)') 'INIT GLEN - FILENAME MISSING'
303 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
304 & SQUEEZE_RIGHT , 1)
305 ENDIF
306
307 ELSE IF (STREAMICEGlenConstConfig.EQ.'UNIFORM' ) THEN
308
309 DO bj = myByLo(myThid), myByHi(myThid)
310 DO bi = myBxLo(myThid), myBxHi(myThid)
311 DO j=1,sNy
312 DO i=1,sNx
313 B_glen(i,j,bi,bj) = B_glen_isothermal
314 ENDDO
315 ENDDO
316 ENDDO
317 ENDDO
318
319 ELSE
320
321 WRITE(msgBuf,'(A)') 'INIT GLEN CONSTANT - NOT IMPLENTED'
322 CALL PRINT_ERROR( msgBuf, myThid)
323 STOP 'ABNORMAL END: S/R STREAMICE_INIT_VAR'
324 ENDIF
325
326 ! finish initialize glen constant
327
328 ! initialize basal traction
329
330 IF ( STREAMICEbasalTracConfig.EQ.'FILE' ) THEN
331
332 IF ( STREAMICEbasalTracFile .NE. ' ' ) THEN
333 _BARRIER
334
335 CALL READ_FLD_XY_RL( STREAMICEbasalTracFile, ' ',
336 & C_basal_friction, 0, myThid )
337
338 ELSE
339 WRITE(msgBuf,'(A)') 'INIT C_BASAL - FILENAME MISSING'
340 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
341 & SQUEEZE_RIGHT , 1)
342 ENDIF
343
344 ELSE IF (STREAMICEbasalTracConfig.EQ.'UNIFORM' ) THEN
345
346 DO bj = myByLo(myThid), myByHi(myThid)
347 DO bi = myBxLo(myThid), myBxHi(myThid)
348 DO j=1,sNy
349 DO i=1,sNx
350 C_basal_friction(i,j,bi,bj) = C_basal_fric_const
351 ENDDO
352 ENDDO
353 ENDDO
354 ENDDO
355
356 ELSE IF (STREAMICEbasalTracConfig.EQ.'2DPERIODIC' ) THEN
357
358 lenx = sNx*nSx*nPx*delX(1)
359 leny = sNy*nSy*nPy*delY(1)
360 DO bj = myByLo(myThid), myByHi(myThid)
361 DO bi = myBxLo(myThid), myBxHi(myThid)
362 DO j=1,sNy
363 DO i=1,sNx
364 x = xC(i,j,bi,bj)
365 y = yC(i,j,bi,bj)
366 C_basal_friction(i,j,bi,bj) =
367 & sqrt(C_basal_fric_const**2*
368 & (1+sin(2*streamice_kx_b_init*PI*x/lenx)*
369 & sin(2*streamice_ky_b_init*PI*y/leny)))
370 ENDDO
371 ENDDO
372 ENDDO
373 ENDDO
374
375 ELSE IF (STREAMICEbasalTracConfig.EQ.'1DPERIODIC' ) THEN
376
377 lenx = sNx*nSx*nPx*delX(1)
378 DO bj = myByLo(myThid), myByHi(myThid)
379 DO bi = myBxLo(myThid), myBxHi(myThid)
380 DO j=1,sNy
381 DO i=1,sNx
382 x = xC(i,j,bi,bj)
383 y = yC(i,j,bi,bj)
384 C_basal_friction(i,j,bi,bj) =
385 & sqrt(C_basal_fric_const**2*(1+
386 & sin(2*streamice_kx_b_init*PI*x/lenx)))
387 ENDDO
388 ENDDO
389 ENDDO
390 ENDDO
391
392 ELSE
393
394 WRITE(msgBuf,'(A)') 'INIT TRAC - NOT IMPLENTED'
395 CALL PRINT_ERROR( msgBuf, myThid)
396 STOP 'ABNORMAL END: S/R STREAMICE_INIT_VAR'
397 ENDIF
398
399 ! finish initialize basal trac
400
401 #ifdef ALLOW_STREAMICE_2DTRACER
402
403 IF ( STREAMICETRAC2DINITFILE .NE. ' ' ) THEN
404 _BARRIER
405
406 CALL READ_FLD_XY_RL( STREAMICETRAC2dInitFile, ' ',
407 & trac2d, 0, myThid )
408
409 ELSE
410 WRITE(msgBuf,'(A)') 'TRAC2dInit - NO FILE SPECIFIED'
411 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
412 & SQUEEZE_RIGHT , 1)
413 DO bj = myByLo(myThid), myByHi(myThid)
414 DO bi = myBxLo(myThid), myBxHi(myThid)
415 DO j=1,sNy
416 DO i=1,sNx
417 trac2d(i,j,bi,bj) = 0.0
418 ENDDO
419 ENDDO
420 ENDDO
421 ENDDO
422
423 ENDIF
424
425 _EXCH_XY_RL (trac2d, myThid)
426
427
428 #endif /*STREAMICE_ALLOW_2DTRACER*/
429
430 CALL STREAMICE_UPD_FFRAC_UNCOUPLED ( myThid )
431
432 _EXCH_XY_RL(H_streamice, myThid )
433 _EXCH_XY_RS(STREAMICE_hmask, myThid )
434 _EXCH_XY_RL(area_shelf_streamice, myThid )
435 _EXCH_XY_RL(C_basal_friction, myThid )
436 _EXCH_XY_RL(B_glen, myThid )
437 #ifdef USE_ALT_RLOW
438 _EXCH_XY_RL(R_low_si, myThid )
439 #endif
440
441 !#ifdef STREAMICE_HYBRID_STRESS
442
443 ! CALL STREAMICE_VISC_BETA (myThid)
444
445 ! DNG THIS CALL IS TO INITIALISE VISCOSITY
446 ! TO AVOID POSSIBLE ADJOINT INSTABILITIES
447 ! IT IS WRITTEN OVER IN FIRST TIMESTEP
448
449 #ifdef ALLOW_AUTODIFF
450
451 CALL STREAMICE_UPD_FFRAC_UNCOUPLED ( myThid )
452 CALL STREAMICE_VELMASK_UPD (myThid)
453 CALL STREAMICE_VEL_SOLVE( myThid )
454
455 #endif
456
457 !#endif
458
459 CALL WRITE_FLD_XY_RL ( "C_basal_fric", "",
460 & C_basal_friction, 0, myThid )
461 CALL WRITE_FLD_XY_RL ( "B_glen_sqrt", "",
462 & B_glen, 0, myThid )
463 CALL WRITE_FLD_XY_RL ( "H_streamice", "init",
464 & H_streamIce, 0, myThid )
465 #ifdef ALLOW_STREAMICE_2DTRACER
466 CALL WRITE_FLD_XY_RL ( "2DTracer", "init",
467 & trac2d, 0, myThid )
468 #endif
469 CALL WRITE_FLD_XY_RL ( "area_shelf_streamice", "init",
470 & area_shelf_streamice, 0, myThid )
471 CALL WRITE_FLD_XY_RL ( "STREAMICE_hmask", "init",
472 & STREAMICE_hmask, 0, myThid )
473 #ifdef ALLOW_CTRL
474 CALL ACTIVE_WRITE_GEN_RS( 'maskCtrlst', STREAMICE_ctrl_mask,
475 & 'XY', Nr, 1, .FALSE., 0, mythid, dummyRS )
476 #endif
477 ! call active_write_xyz( 'maskCtrlS', STREAMICE_ctrl_mask, 1, 0,
478 ! & mythid, dummy)
479 ! CALL STREAMICE_VELMASK_UPD (myThid)
480 ! CALL STREAMICE_UPD_FFRAC_UNCOUPLED ( myThid )
481 ! CALL STREAMICE_VEL_SOLVE( myThid )
482
483 CALL WRITE_FLD_XY_RL ( "U_init", "",
484 & C_basal_friction, 0, myThid )
485 CALL WRITE_FLD_XY_RL ( "V_init", "",
486 & V_streamice, 0, myThid )
487 #ifdef USE_ALT_RLOW
488 CALL WRITE_FLD_XY_RL ( "R_low_si", "init",
489 & R_low_si, 0, myThid )
490 #endif
491
492 ! CALL WRITE_FULLARRAY_RL ("H",H_streamice,1,0,0,1,0,myThid)
493 ! CALL WRITE_FULLARRAY_RL ("hmask",STREAMICE_hmask,1,0,0,1,0,myThid)
494 ! CALL WRITE_FULLARRAY_RL ("umask",STREAMICE_umask,1,0,0,1,0,myThid)
495
496
497
498 #endif /* ALLOW_STREAMICE */
499
500 RETURN
501 END
502

  ViewVC Help
Powered by ViewVC 1.1.22