/[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.6 - (show annotations) (download)
Wed Mar 26 17:40:11 2014 UTC (10 years, 2 months ago) by dgoldberg
Branch: MAIN
CVS Tags: checkpoint64v
Changes since 1.5: +10 -1 lines
streamice pickup; impose calving boundary beyond init. front

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

  ViewVC Help
Powered by ViewVC 1.1.22