/[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.4 - (show annotations) (download)
Fri Jun 21 22:00:03 2013 UTC (10 years, 11 months ago) by jmc
Branch: MAIN
Changes since 1.3: +39 -40 lines
fix other RS/RL mismatch

1 C $Header: /u/gcmpack/MITgcm/pkg/streamice/streamice_init_varia.F,v 1.3 2013/06/21 21:15:36 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 ! finish initialize thickness
291
292 ! initialize glen constant
293
294 IF ( STREAMICEGlenConstConfig.EQ.'FILE' ) THEN
295
296 IF ( STREAMICEGlenConstFile .NE. ' ' ) THEN
297 _BARRIER
298
299 CALL READ_FLD_XY_RL( STREAMICEGlenConstFile, ' ',
300 & B_glen, 0, myThid )
301
302 ELSE
303 WRITE(msgBuf,'(A)') 'INIT GLEN - FILENAME MISSING'
304 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
305 & SQUEEZE_RIGHT , 1)
306 ENDIF
307
308 ELSE IF (STREAMICEGlenConstConfig.EQ.'UNIFORM' ) THEN
309
310 DO bj = myByLo(myThid), myByHi(myThid)
311 DO bi = myBxLo(myThid), myBxHi(myThid)
312 DO j=1,sNy
313 DO i=1,sNx
314 B_glen(i,j,bi,bj) = B_glen_isothermal
315 ENDDO
316 ENDDO
317 ENDDO
318 ENDDO
319
320 ELSE
321
322 WRITE(msgBuf,'(A)') 'INIT GLEN CONSTANT - NOT IMPLENTED'
323 CALL PRINT_ERROR( msgBuf, myThid)
324 STOP 'ABNORMAL END: S/R STREAMICE_INIT_VAR'
325 ENDIF
326
327 ! finish initialize glen constant
328
329 ! initialize basal traction
330
331 IF ( STREAMICEbasalTracConfig.EQ.'FILE' ) THEN
332
333 IF ( STREAMICEbasalTracFile .NE. ' ' ) THEN
334 _BARRIER
335
336 CALL READ_FLD_XY_RL( STREAMICEbasalTracFile, ' ',
337 & C_basal_friction, 0, myThid )
338
339 ELSE
340 WRITE(msgBuf,'(A)') 'INIT C_BASAL - FILENAME MISSING'
341 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
342 & SQUEEZE_RIGHT , 1)
343 ENDIF
344
345 ELSE IF (STREAMICEbasalTracConfig.EQ.'UNIFORM' ) THEN
346
347 DO bj = myByLo(myThid), myByHi(myThid)
348 DO bi = myBxLo(myThid), myBxHi(myThid)
349 DO j=1,sNy
350 DO i=1,sNx
351 C_basal_friction(i,j,bi,bj) = C_basal_fric_const
352 ENDDO
353 ENDDO
354 ENDDO
355 ENDDO
356
357 ELSE IF (STREAMICEbasalTracConfig.EQ.'2DPERIODIC' ) THEN
358
359 lenx = sNx*nSx*nPx*delX(1)
360 leny = sNy*nSy*nPy*delY(1)
361 DO bj = myByLo(myThid), myByHi(myThid)
362 DO bi = myBxLo(myThid), myBxHi(myThid)
363 DO j=1,sNy
364 DO i=1,sNx
365 x = xC(i,j,bi,bj)
366 y = yC(i,j,bi,bj)
367 C_basal_friction(i,j,bi,bj) =
368 & sqrt(C_basal_fric_const**2*
369 & (1+sin(2*streamice_kx_b_init*PI*x/lenx)*
370 & sin(2*streamice_ky_b_init*PI*y/leny)))
371 ENDDO
372 ENDDO
373 ENDDO
374 ENDDO
375
376 ELSE IF (STREAMICEbasalTracConfig.EQ.'1DPERIODIC' ) THEN
377
378 lenx = sNx*nSx*nPx*delX(1)
379 DO bj = myByLo(myThid), myByHi(myThid)
380 DO bi = myBxLo(myThid), myBxHi(myThid)
381 DO j=1,sNy
382 DO i=1,sNx
383 x = xC(i,j,bi,bj)
384 y = yC(i,j,bi,bj)
385 C_basal_friction(i,j,bi,bj) =
386 & sqrt(C_basal_fric_const**2*(1+
387 & sin(2*streamice_kx_b_init*PI*x/lenx)))
388 ENDDO
389 ENDDO
390 ENDDO
391 ENDDO
392
393 ELSE
394
395 WRITE(msgBuf,'(A)') 'INIT TRAC - NOT IMPLENTED'
396 CALL PRINT_ERROR( msgBuf, myThid)
397 STOP 'ABNORMAL END: S/R STREAMICE_INIT_VAR'
398 ENDIF
399
400 ! finish initialize basal trac
401
402 #ifdef ALLOW_STREAMICE_2DTRACER
403
404 IF ( STREAMICETRAC2DINITFILE .NE. ' ' ) THEN
405 _BARRIER
406
407 CALL READ_FLD_XY_RL( STREAMICETRAC2dInitFile, ' ',
408 & trac2d, 0, myThid )
409
410 ELSE
411 WRITE(msgBuf,'(A)') 'TRAC2dInit - NO FILE SPECIFIED'
412 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
413 & SQUEEZE_RIGHT , 1)
414 DO bj = myByLo(myThid), myByHi(myThid)
415 DO bi = myBxLo(myThid), myBxHi(myThid)
416 DO j=1,sNy
417 DO i=1,sNx
418 trac2d(i,j,bi,bj) = 0.0
419 ENDDO
420 ENDDO
421 ENDDO
422 ENDDO
423
424 ENDIF
425
426 _EXCH_XY_RL (trac2d, myThid)
427
428
429 #endif /*STREAMICE_ALLOW_2DTRACER*/
430
431 CALL STREAMICE_UPD_FFRAC_UNCOUPLED ( myThid )
432
433 _EXCH_XY_RL(H_streamice, myThid )
434 _EXCH_XY_RS(STREAMICE_hmask, myThid )
435 _EXCH_XY_RL(area_shelf_streamice, myThid )
436 _EXCH_XY_RL(C_basal_friction, myThid )
437 _EXCH_XY_RL(B_glen, myThid )
438 #ifdef USE_ALT_RLOW
439 _EXCH_XY_RL(R_low_si, myThid )
440 #endif
441
442 !#ifdef STREAMICE_HYBRID_STRESS
443
444 ! CALL STREAMICE_VISC_BETA (myThid)
445
446 ! DNG THIS CALL IS TO INITIALISE VISCOSITY
447 ! TO AVOID POSSIBLE ADJOINT INSTABILITIES
448 ! IT IS WRITTEN OVER IN FIRST TIMESTEP
449
450 #ifdef ALLOW_AUTODIFF
451
452 CALL STREAMICE_UPD_FFRAC_UNCOUPLED ( myThid )
453 CALL STREAMICE_VELMASK_UPD (myThid)
454 CALL STREAMICE_VEL_SOLVE( myThid )
455
456 #endif
457
458 !#endif
459
460 CALL WRITE_FLD_XY_RL ( "C_basal_fric", "",
461 & C_basal_friction, 0, myThid )
462 CALL WRITE_FLD_XY_RL ( "B_glen_sqrt", "",
463 & B_glen, 0, myThid )
464 CALL WRITE_FLD_XY_RL ( "H_streamice", "init",
465 & H_streamIce, 0, myThid )
466 #ifdef ALLOW_STREAMICE_2DTRACER
467 CALL WRITE_FLD_XY_RL ( "2DTracer", "init",
468 & trac2d, 0, myThid )
469 #endif
470 CALL WRITE_FLD_XY_RL ( "area_shelf_streamice", "init",
471 & area_shelf_streamice, 0, myThid )
472 CALL WRITE_FLD_XY_RS ( "STREAMICE_hmask", "init",
473 & STREAMICE_hmask, 0, myThid )
474 #ifdef ALLOW_CTRL
475 CALL ACTIVE_WRITE_GEN_RS( 'maskCtrlst', STREAMICE_ctrl_mask,
476 & 'XY', Nr, 1, .FALSE., 0, myThid, dummyRS )
477 #endif
478 ! call active_write_xyz( 'maskCtrlS', STREAMICE_ctrl_mask, 1, 0,
479 ! & myThid, dummy)
480 ! CALL STREAMICE_VELMASK_UPD (myThid)
481 ! CALL STREAMICE_UPD_FFRAC_UNCOUPLED ( myThid )
482 ! CALL STREAMICE_VEL_SOLVE( myThid )
483
484 CALL WRITE_FLD_XY_RL ( "U_init", "",
485 & C_basal_friction, 0, myThid )
486 CALL WRITE_FLD_XY_RL ( "V_init", "",
487 & V_streamice, 0, myThid )
488 #ifdef USE_ALT_RLOW
489 CALL WRITE_FLD_XY_RL ( "R_low_si", "init",
490 & R_low_si, 0, myThid )
491 #endif
492
493 ! CALL WRITE_FULLARRAY_RL ("H",H_streamice,1,0,0,1,0,myThid)
494 ! CALL WRITE_FULLARRAY_RS ("hmask",STREAMICE_hmask,1,0,0,1,0,myThid)
495 ! CALL WRITE_FULLARRAY_RL ("umask",STREAMICE_umask,1,0,0,1,0,myThid)
496
497 #endif /* ALLOW_STREAMICE */
498
499 RETURN
500 END
501

  ViewVC Help
Powered by ViewVC 1.1.22