/[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.13 - (show annotations) (download)
Thu May 23 22:12:33 2013 UTC (12 years, 2 months ago) by dgoldberg
Branch: MAIN
Changes since 1.12: +51 -1 lines
specify melt rate field

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

  ViewVC Help
Powered by ViewVC 1.1.22