/[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.5 - (show annotations) (download)
Thu Sep 20 02:04:45 2012 UTC (12 years, 10 months ago) by heimbach
Branch: MAIN
Changes since 1.4: +4 -2 lines
Enable working version for #define STREAMICE_HYBRID_STRESS

1 C $Header: /u/gcmpack/MITgcm_contrib/dgoldberg/streamice/streamice_init_varia.F,v 1.4 2012/09/18 17:06:48 dgoldberg Exp $
2 C $Name: $
3
4 #include "STREAMICE_OPTIONS.h"
5
6 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
7
8 CBOP
9 SUBROUTINE STREAMICE_INIT_VARIA( myThid )
10 C /============================================================\
11 C | SUBROUTINE STREAMICE_INIT_VARIA |
12 C | o Routine to initialize STREAMICE variables. |
13 C |============================================================|
14 C | Initialize STREAMICE parameters and variables. |
15 C \============================================================/
16 IMPLICIT NONE
17
18 C === Global variables ===
19 #include "SIZE.h"
20 #include "GRID.h"
21 #include "SET_GRID.h"
22 #include "EEPARAMS.h"
23 #include "PARAMS.h"
24 #include "STREAMICE.h"
25 #include "STREAMICE_CG.h"
26 #include "STREAMICE_ADV.h"
27
28 C === Routine arguments ===
29 C myThid - Number of this instance of STREAMICE_INIT_VARIA
30 INTEGER myThid
31 CEndOfInterface
32
33 #ifdef ALLOW_STREAMICE
34 C === Local variables ===
35 C I,J,bi,bj - Loop counters
36 INTEGER i, j, k, bi, bj, Gi, Gj
37 INTEGER col_y, col_x
38 _RL slope_pos, c1, x, y, lenx, leny
39 CHARACTER*(MAX_LEN_MBUF) msgBuf
40 CEOP
41
42 C ZERO OUT FLOATING POINT ARRAYS
43
44 DO bj = myByLo(myThid), myByHi(myThid)
45 DO bi = myBxLo(myThid), myBxHi(myThid)
46 DO j=1-Oly,sNy+Oly
47 DO i=1-Olx,sNx+Olx
48 H_streamIce(i,j,bi,bj) = 0. _d 0
49 U_streamice(i,j,bi,bj) = 0. _d 0
50 V_streamice(i,j,bi,bj) = 0. _d 0
51 visc_streamice(i,j,bi,bj) = 0. _d 0
52 tau_beta_eff_streamice(i,j,bi,bj) = 0. _d 0
53 float_frac_streamice(i,j,bi,bj) = 0. _d 0
54 base_el_streamice(i,j,bi,bj) = 0. _d 0
55 surf_el_streamice(i,j,bi,bj) = 0. _d 0
56 area_shelf_streamice(i,j,bi,bj) = 0. _d 0
57 mass_ice_streamice(i,j,bi,bj) = 0. _d 0
58 BDOT_streamice(i,j,bi,bj) = 0. _d 0
59 ! C_basal_friction(i,j,bi,bj) = C_basal_fric_const
60 A_glen(i,j,bi,bj) = A_glen_isothermal
61 #ifdef ALLOW_AUTODIFF_TAMC
62 ru_old_si(i,j,bi,bj) = 0. _d 0
63 rv_old_si(i,j,bi,bj) = 0. _d 0
64 zu_old_si(i,j,bi,bj) = 0. _d 0
65 zv_old_si(i,j,bi,bj) = 0. _d 0
66 h_after_uflux_SI(i,j,bi,bj) = 0. _d 0
67 #endif
68 #ifdef STREAMICE_HYBRID_STRESS
69 do k=1,Nr
70 visc_streamice_full(i,j,k,bi,bj) =
71 & eps_glen_min**((1-n_glen)/n_glen)
72 enddo
73 streamice_taubx (i,j,bi,bj) = 0. _d 0
74 streamice_tauby (i,j,bi,bj) = 0. _d 0
75 #endif
76 ENDDO
77 ENDDO
78 ENDDO
79 ENDDO
80
81 DO j = 1-oly, sNy+oly
82 DO i = 1-olx, sNx+olx
83 DO bj = myByLo(myThid), myByHi(myThid)
84 DO bi = myBxLo(myThid), myBxHi(myThid)
85 cc DO k=1,4
86 DO col_x=-1,1
87 DO col_y=-1,1
88 streamice_cg_A1(i,j,bi,bj,col_x,col_y)=0.0
89 streamice_cg_A2(i,j,bi,bj,col_x,col_y)=0.0
90 streamice_cg_A3(i,j,bi,bj,col_x,col_y)=0.0
91 streamice_cg_A4(i,j,bi,bj,col_x,col_y)=0.0
92 ENDDO
93 ENDDO
94 cc ENDDO
95 ENDDO
96 ENDDO
97 ENDDO
98 ENDDO
99
100 C INIT. INTEGER ARRAYS
101
102 DO bj = myByLo(myThid), myByHi(myThid)
103 DO bi = myBxLo(myThid), myBxHi(myThid)
104 DO j=1-Oly,sNy+Oly
105 DO i=1-Olx,sNx+Olx
106 STREAMICE_hmask(i,j,bi,bj) = -1.0
107 STREAMICE_umask(i,j,bi,bj) = 0.0
108 STREAMICE_vmask(i,j,bi,bj) = 0.0
109 STREAMICE_ufacemask(i,j,bi,bj) = 0.0
110 STREAMICE_vfacemask(i,j,bi,bj) = 0.0
111 STREAMICE_float_cond(i,j,bi,bj) = 0.0
112 ENDDO
113 ENDDO
114 ENDDO
115 ENDDO
116
117 ! initialize thickness
118
119 IF ( STREAMICEthickInit.EQ.'PARAM' ) THEN
120
121 WRITE(msgBuf,'(A)') 'initializing analytic thickness'
122 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
123 & SQUEEZE_RIGHT , 1)
124
125 slope_pos = shelf_edge_pos - shelf_flat_width
126 c1 = 0.0
127 IF (shelf_slope_scale .GT. 0.0) THEN
128 c1 = 1.0 / shelf_slope_scale
129 ENDIF
130
131 DO bj = myByLo(myThid), myByHi(myThid)
132 DO bi = myBxLo(myThid), myBxHi(myThid)
133 DO j=1,sNy
134 DO i=1,sNx
135 Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
136 Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
137
138 IF ((Gi.lt.Nx).and.(Gj.lt.Ny)) THEN
139
140 C IF (flow_dir .EQ. 2.0) THEN
141 IF (.TRUE.) THEN
142 IF (xC(i-1,j,bi,bj).GE.shelf_edge_pos) THEN
143 area_shelf_streamice(i,j,bi,bj) = 0. _d 0
144 STREAMICE_hmask(i,j,bi,bj) = 0. _d 0
145 ELSE
146
147 IF (xC(i,j,bi,bj).GT.slope_pos) THEN
148 H_streamice (i,j,bi,bj) = shelf_min_draft
149 ELSE
150 H_streamice (i,j,bi,bj) = (shelf_min_draft +
151 & (shelf_max_draft - shelf_min_draft) *
152 & min (1.0, (c1*(slope_pos-xC(i,j,bi,bj)))**2))
153 ENDIF
154
155 IF (xC(i,j,bi,bj).GT.shelf_edge_pos) THEN
156 area_shelf_streamice(i,j,bi,bj) = rA(i,j,bi,bj) *
157 & (shelf_edge_pos-xG(i,j,bi,bj)) /
158 & (xG(i+1,j,bi,bj)-xG(i,j,bi,bj))
159 IF (area_shelf_streamice(i,j,bi,bj).gt. 0. _d 0) THEN
160 STREAMICE_hmask(i,j,bi,bj) = 2.0
161 ELSE
162 STREAMICE_hmask(i,j,bi,bj) = 0.0
163 H_streamice(i,j,bi,bj) = 0.0
164 ENDIF
165 ELSE
166 area_shelf_streamice(i,j,bi,bj) = rA(i,j,bi,bj)
167 STREAMICE_hmask(i,j,bi,bj) = 1.0
168 ENDIF
169
170
171 ENDIF
172 ENDIF
173 ENDIF
174 ENDDO
175 ENDDO
176 ENDDO
177 ENDDO
178
179 ELSE IF ( STREAMICEthickInit.EQ.'FILE' ) THEN
180
181 IF ( STREAMICEthickFile .NE. ' ' ) THEN
182 _BARRIER
183 C The 0 is the "iteration" argument. The ' ' is an empty suffix
184 CALL READ_FLD_XY_RS( STREAMICEthickFile, ' ', H_streamice,
185 & 0, myThid )
186 DO bj = myByLo(myThid), myByHi(myThid)
187 DO bi = myBxLo(myThid), myBxHi(myThid)
188 DO j=1,sNy
189 DO i=1,sNx
190 Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
191 Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
192 IF ((Gi.lt.Nx.OR.STREAMICE_EW_periodic).and.
193 & (Gj.lt.Ny.OR.STREAMICE_NS_periodic)) THEN
194 IF (H_streamice(i,j,bi,bj).GT.0. _d 0) THEN
195 area_shelf_streamice(i,j,bi,bj) = rA(i,j,bi,bj)
196 STREAMICE_hmask(i,j,bi,bj) = 1.0
197 ELSE
198 area_shelf_streamice(i,j,bi,bj) = 0. _d 0
199 STREAMICE_hmask(i,j,bi,bj) = 0. _d 0
200 ENDIF
201 ENDIF
202 ENDDO
203 ENDDO
204 ENDDO
205 ENDDO
206 ELSE
207 WRITE(msgBuf,'(A)') 'INIT THICKNESS - FILENAME MISSING'
208 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
209 & SQUEEZE_RIGHT , 1)
210 ENDIF
211
212 ELSE
213
214 WRITE(msgBuf,'(A)') 'INIT THICKNESS - NOT IMPLENTED'
215 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
216 & SQUEEZE_RIGHT , 1)
217 ENDIF
218
219 ! finish initialize thickness
220
221 ! initialize basal traction
222
223 IF ( STREAMICEbasalTracConfig.EQ.'FILE' ) THEN
224
225 IF ( STREAMICEbasalTracFile .NE. ' ' ) THEN
226 _BARRIER
227 C The 0 is the "iteration" argument. The ' ' is an empty suffix
228 CALL READ_FLD_XY_RS( STREAMICEbasalTracFile, ' ',
229 & C_basal_friction, 0, myThid )
230
231 ELSE
232 WRITE(msgBuf,'(A)') 'INIT THICKNESS - FILENAME MISSING'
233 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
234 & SQUEEZE_RIGHT , 1)
235 ENDIF
236
237 ELSE IF (STREAMICEbasalTracConfig.EQ.'UNIFORM' ) THEN
238
239 DO bj = myByLo(myThid), myByHi(myThid)
240 DO bi = myBxLo(myThid), myBxHi(myThid)
241 DO j=1,sNy
242 DO i=1,sNx
243 C_basal_friction(i,j,bi,bj) = C_basal_fric_const
244 ENDDO
245 ENDDO
246 ENDDO
247 ENDDO
248
249 ELSE IF (STREAMICEbasalTracConfig.EQ.'2DPERIODIC' ) THEN
250
251 lenx = sNx*nSx*nPx*delX(1)
252 leny = sNy*nSy*nPy*delY(1)
253 print *, 'lenx', lenx
254 print *, 'leny', leny
255 DO bj = myByLo(myThid), myByHi(myThid)
256 DO bi = myBxLo(myThid), myBxHi(myThid)
257 DO j=1,sNy
258 DO i=1,sNx
259 x = xC(i,j,bi,bj)
260 y = yC(i,j,bi,bj)
261 C_basal_friction(i,j,bi,bj) =
262 & sqrt(C_basal_fric_const**2*
263 & (1+sin(2*streamice_kx_b_init*PI*x/lenx)*
264 & sin(2*streamice_ky_b_init*PI*y/leny)))
265 ENDDO
266 ENDDO
267 ENDDO
268 ENDDO
269
270 ELSE IF (STREAMICEbasalTracConfig.EQ.'1DPERIODIC' ) THEN
271
272 lenx = sNx*nSx*nPx*delX(1)
273 print *, 'lenx', lenx
274 DO bj = myByLo(myThid), myByHi(myThid)
275 DO bi = myBxLo(myThid), myBxHi(myThid)
276 DO j=1,sNy
277 DO i=1,sNx
278 x = xC(i,j,bi,bj)
279 y = yC(i,j,bi,bj)
280 C_basal_friction(i,j,bi,bj) =
281 & sqrt(C_basal_fric_const**2*(1+
282 & sin(2*streamice_kx_b_init*PI*x/lenx)))
283 ENDDO
284 ENDDO
285 ENDDO
286 ENDDO
287
288 ELSE
289
290 WRITE(msgBuf,'(A)') 'INIT THICKNESS - NOT IMPLENTED'
291 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
292 & SQUEEZE_RIGHT , 1)
293 ENDIF
294
295 ! finish initialize basal traction
296
297 CALL STREAMICE_UPD_FFRAC_UNCOUPLED ( myThid )
298
299 _EXCH_XY_RL(H_streamice, myThid )
300 _EXCH_XY_RL(STREAMICE_hmask, myThid )
301 _EXCH_XY_RL(area_shelf_streamice, myThid )
302 _EXCH_XY_RL(C_basal_friction, myThid )
303
304
305 #ifdef STREAMICE_HYBRID_STRESS
306
307 CALL STREAMICE_VISC_BETA (myThid)
308
309 #endif
310
311 CALL WRITE_FLD_XY_RL ( "H_streamice", "init",
312 & H_streamIce, 0, myThid )
313 CALL WRITE_FLD_XY_RL ( "area_shelf_streamice", "init",
314 & area_shelf_streamice, 0, myThid )
315 CALL WRITE_FLD_XY_RL ( "STREAMICE_hmask", "init",
316 & STREAMICE_hmask, 0, myThid )
317
318 ! CALL STREAMICE_VELMASK_UPD (myThid)
319 ! CALL STREAMICE_UPD_FFRAC_UNCOUPLED ( myThid )
320 ! CALL STREAMICE_VEL_SOLVE( myThid )
321
322 CALL WRITE_FLD_XY_RL ( "U_init", "",
323 & C_basal_friction, 0, myThid )
324 CALL WRITE_FLD_XY_RL ( "V_init", "",
325 & V_streamice, 0, myThid )
326
327 ! CALL WRITE_FULLARRAY_RL ("H",H_streamice,1,0,0,1,0,myThid)
328 ! CALL WRITE_FULLARRAY_RL ("hmask",STREAMICE_hmask,1,0,0,1,0,myThid)
329 ! CALL WRITE_FULLARRAY_RL ("umask",STREAMICE_umask,1,0,0,1,0,myThid)
330
331 #endif /* ALLOW_STREAMICE */
332
333 RETURN
334 END
335

  ViewVC Help
Powered by ViewVC 1.1.22