/[MITgcm]/MITgcm_contrib/dgoldberg/streamice/streamice_init_varia.F
ViewVC logotype

Annotation of /MITgcm_contrib/dgoldberg/streamice/streamice_init_varia.F

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


Revision 1.9 - (hide annotations) (download)
Wed Nov 21 01:04:21 2012 UTC (12 years, 8 months ago) by dgoldberg
Branch: MAIN
Changes since 1.8: +29 -6 lines
call to solve velocity (and initialize ice viscosity) if AD is turned on; velocities overwritten in first timestep

1 dgoldberg 1.9 C $Header: /u/gcmpack/MITgcm_contrib/dgoldberg/streamice/streamice_init_varia.F,v 1.8 2012/10/10 15:03:10 dgoldberg Exp $
2 heimbach 1.1 C $Name: $
3 dgoldberg 1.9 #ifdef ALLOW_COST
4     # include "COST_OPTIONS.h"
5     #endif
6 heimbach 1.1 #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 dgoldberg 1.4 #include "SET_GRID.h"
24 heimbach 1.1 #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 dgoldberg 1.4 _RL slope_pos, c1, x, y, lenx, leny
41 heimbach 1.1 CHARACTER*(MAX_LEN_MBUF) msgBuf
42 dgoldberg 1.9 _RS dummyRS
43    
44 heimbach 1.1 CEOP
45    
46     C ZERO OUT FLOATING POINT ARRAYS
47    
48     DO bj = myByLo(myThid), myByHi(myThid)
49     DO bi = myBxLo(myThid), myBxHi(myThid)
50     DO j=1-Oly,sNy+Oly
51     DO i=1-Olx,sNx+Olx
52     H_streamIce(i,j,bi,bj) = 0. _d 0
53     U_streamice(i,j,bi,bj) = 0. _d 0
54     V_streamice(i,j,bi,bj) = 0. _d 0
55     visc_streamice(i,j,bi,bj) = 0. _d 0
56     tau_beta_eff_streamice(i,j,bi,bj) = 0. _d 0
57     float_frac_streamice(i,j,bi,bj) = 0. _d 0
58     base_el_streamice(i,j,bi,bj) = 0. _d 0
59     surf_el_streamice(i,j,bi,bj) = 0. _d 0
60     area_shelf_streamice(i,j,bi,bj) = 0. _d 0
61     mass_ice_streamice(i,j,bi,bj) = 0. _d 0
62 dgoldberg 1.3 BDOT_streamice(i,j,bi,bj) = 0. _d 0
63 dgoldberg 1.4 ! C_basal_friction(i,j,bi,bj) = C_basal_fric_const
64 dgoldberg 1.3 A_glen(i,j,bi,bj) = A_glen_isothermal
65 dgoldberg 1.7 H_streamice_prev(i,j,bi,bj) = 0. _d 0
66 heimbach 1.1 #ifdef ALLOW_AUTODIFF_TAMC
67     ru_old_si(i,j,bi,bj) = 0. _d 0
68     rv_old_si(i,j,bi,bj) = 0. _d 0
69     zu_old_si(i,j,bi,bj) = 0. _d 0
70     zv_old_si(i,j,bi,bj) = 0. _d 0
71     h_after_uflux_SI(i,j,bi,bj) = 0. _d 0
72     #endif
73 dgoldberg 1.6 #ifdef USE_ALT_RLOW
74     R_low_si(i,j,bi,bj) = 0. _d 0
75     #endif
76    
77 dgoldberg 1.4 #ifdef STREAMICE_HYBRID_STRESS
78     do k=1,Nr
79     visc_streamice_full(i,j,k,bi,bj) =
80     & eps_glen_min**((1-n_glen)/n_glen)
81     enddo
82 heimbach 1.5 streamice_taubx (i,j,bi,bj) = 0. _d 0
83     streamice_tauby (i,j,bi,bj) = 0. _d 0
84 dgoldberg 1.4 #endif
85 heimbach 1.1 ENDDO
86     ENDDO
87 dgoldberg 1.6
88     #ifdef ALLOW_COST_TEST
89     cost_func1_streamice (bi,bj) = 0.0
90     #endif
91    
92 heimbach 1.1 ENDDO
93     ENDDO
94    
95     DO j = 1-oly, sNy+oly
96     DO i = 1-olx, sNx+olx
97     DO bj = myByLo(myThid), myByHi(myThid)
98     DO bi = myBxLo(myThid), myBxHi(myThid)
99     cc DO k=1,4
100     DO col_x=-1,1
101     DO col_y=-1,1
102     streamice_cg_A1(i,j,bi,bj,col_x,col_y)=0.0
103     streamice_cg_A2(i,j,bi,bj,col_x,col_y)=0.0
104     streamice_cg_A3(i,j,bi,bj,col_x,col_y)=0.0
105     streamice_cg_A4(i,j,bi,bj,col_x,col_y)=0.0
106     ENDDO
107     ENDDO
108     cc ENDDO
109     ENDDO
110     ENDDO
111     ENDDO
112     ENDDO
113    
114     C INIT. INTEGER ARRAYS
115    
116     DO bj = myByLo(myThid), myByHi(myThid)
117     DO bi = myBxLo(myThid), myBxHi(myThid)
118     DO j=1-Oly,sNy+Oly
119     DO i=1-Olx,sNx+Olx
120     STREAMICE_hmask(i,j,bi,bj) = -1.0
121     STREAMICE_umask(i,j,bi,bj) = 0.0
122     STREAMICE_vmask(i,j,bi,bj) = 0.0
123     STREAMICE_ufacemask(i,j,bi,bj) = 0.0
124     STREAMICE_vfacemask(i,j,bi,bj) = 0.0
125     STREAMICE_float_cond(i,j,bi,bj) = 0.0
126     ENDDO
127     ENDDO
128     ENDDO
129     ENDDO
130    
131 dgoldberg 1.6
132     #ifdef USE_ALT_RLOW
133     ! init alternate array for topog
134     IF ( bathyFile .NE. ' ' ) THEN
135     _BARRIER
136     C The 0 is the "iteration" argument. The ' ' is an empty suffix
137     CALL READ_FLD_XY_RS( bathyFile, '',
138     & R_low_si, 0, myThid )
139    
140     ELSE
141     WRITE(msgBuf,'(A)') 'STREAMICE TOPOG - FILENAME MISSING'
142     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
143     & SQUEEZE_RIGHT , 1)
144     ENDIF
145     #endif
146    
147 dgoldberg 1.4 ! initialize thickness
148 heimbach 1.1
149     IF ( STREAMICEthickInit.EQ.'PARAM' ) THEN
150    
151     WRITE(msgBuf,'(A)') 'initializing analytic thickness'
152     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
153     & SQUEEZE_RIGHT , 1)
154    
155     slope_pos = shelf_edge_pos - shelf_flat_width
156     c1 = 0.0
157     IF (shelf_slope_scale .GT. 0.0) THEN
158     c1 = 1.0 / shelf_slope_scale
159     ENDIF
160    
161     DO bj = myByLo(myThid), myByHi(myThid)
162     DO bi = myBxLo(myThid), myBxHi(myThid)
163     DO j=1,sNy
164     DO i=1,sNx
165     Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
166     Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
167    
168     IF ((Gi.lt.Nx).and.(Gj.lt.Ny)) THEN
169    
170     C IF (flow_dir .EQ. 2.0) THEN
171     IF (.TRUE.) THEN
172     IF (xC(i-1,j,bi,bj).GE.shelf_edge_pos) THEN
173     area_shelf_streamice(i,j,bi,bj) = 0. _d 0
174     STREAMICE_hmask(i,j,bi,bj) = 0. _d 0
175     ELSE
176    
177     IF (xC(i,j,bi,bj).GT.slope_pos) THEN
178     H_streamice (i,j,bi,bj) = shelf_min_draft
179     ELSE
180     H_streamice (i,j,bi,bj) = (shelf_min_draft +
181     & (shelf_max_draft - shelf_min_draft) *
182     & min (1.0, (c1*(slope_pos-xC(i,j,bi,bj)))**2))
183     ENDIF
184    
185     IF (xC(i,j,bi,bj).GT.shelf_edge_pos) THEN
186     area_shelf_streamice(i,j,bi,bj) = rA(i,j,bi,bj) *
187     & (shelf_edge_pos-xG(i,j,bi,bj)) /
188     & (xG(i+1,j,bi,bj)-xG(i,j,bi,bj))
189 heimbach 1.5 IF (area_shelf_streamice(i,j,bi,bj).gt. 0. _d 0) THEN
190 heimbach 1.1 STREAMICE_hmask(i,j,bi,bj) = 2.0
191     ELSE
192     STREAMICE_hmask(i,j,bi,bj) = 0.0
193     H_streamice(i,j,bi,bj) = 0.0
194     ENDIF
195     ELSE
196     area_shelf_streamice(i,j,bi,bj) = rA(i,j,bi,bj)
197     STREAMICE_hmask(i,j,bi,bj) = 1.0
198     ENDIF
199    
200    
201     ENDIF
202     ENDIF
203     ENDIF
204     ENDDO
205     ENDDO
206     ENDDO
207     ENDDO
208    
209     ELSE IF ( STREAMICEthickInit.EQ.'FILE' ) THEN
210    
211     IF ( STREAMICEthickFile .NE. ' ' ) THEN
212     _BARRIER
213     C The 0 is the "iteration" argument. The ' ' is an empty suffix
214     CALL READ_FLD_XY_RS( STREAMICEthickFile, ' ', H_streamice,
215     & 0, myThid )
216     DO bj = myByLo(myThid), myByHi(myThid)
217     DO bi = myBxLo(myThid), myBxHi(myThid)
218     DO j=1,sNy
219     DO i=1,sNx
220     Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
221     Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
222 dgoldberg 1.4 IF ((Gi.lt.Nx.OR.STREAMICE_EW_periodic).and.
223     & (Gj.lt.Ny.OR.STREAMICE_NS_periodic)) THEN
224 heimbach 1.1 IF (H_streamice(i,j,bi,bj).GT.0. _d 0) THEN
225     area_shelf_streamice(i,j,bi,bj) = rA(i,j,bi,bj)
226     STREAMICE_hmask(i,j,bi,bj) = 1.0
227     ELSE
228     area_shelf_streamice(i,j,bi,bj) = 0. _d 0
229     STREAMICE_hmask(i,j,bi,bj) = 0. _d 0
230     ENDIF
231 dgoldberg 1.9 Do k=1,Nr
232     STREAMICE_ctrl_mask(i,j,bi,bj,k) = 1. _d 0
233     enddo
234 heimbach 1.1 ENDIF
235     ENDDO
236     ENDDO
237     ENDDO
238     ENDDO
239     ELSE
240     WRITE(msgBuf,'(A)') 'INIT THICKNESS - FILENAME MISSING'
241     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
242     & SQUEEZE_RIGHT , 1)
243     ENDIF
244    
245     ELSE
246    
247     WRITE(msgBuf,'(A)') 'INIT THICKNESS - NOT IMPLENTED'
248     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
249     & SQUEEZE_RIGHT , 1)
250     ENDIF
251    
252 dgoldberg 1.4 ! finish initialize thickness
253    
254     ! initialize basal traction
255    
256     IF ( STREAMICEbasalTracConfig.EQ.'FILE' ) THEN
257    
258     IF ( STREAMICEbasalTracFile .NE. ' ' ) THEN
259     _BARRIER
260     C The 0 is the "iteration" argument. The ' ' is an empty suffix
261     CALL READ_FLD_XY_RS( STREAMICEbasalTracFile, ' ',
262     & C_basal_friction, 0, myThid )
263    
264     ELSE
265     WRITE(msgBuf,'(A)') 'INIT THICKNESS - FILENAME MISSING'
266     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
267     & SQUEEZE_RIGHT , 1)
268     ENDIF
269    
270     ELSE IF (STREAMICEbasalTracConfig.EQ.'UNIFORM' ) THEN
271    
272     DO bj = myByLo(myThid), myByHi(myThid)
273     DO bi = myBxLo(myThid), myBxHi(myThid)
274     DO j=1,sNy
275     DO i=1,sNx
276     C_basal_friction(i,j,bi,bj) = C_basal_fric_const
277     ENDDO
278     ENDDO
279     ENDDO
280     ENDDO
281    
282     ELSE IF (STREAMICEbasalTracConfig.EQ.'2DPERIODIC' ) THEN
283    
284     lenx = sNx*nSx*nPx*delX(1)
285     leny = sNy*nSy*nPy*delY(1)
286     print *, 'lenx', lenx
287     print *, 'leny', leny
288     DO bj = myByLo(myThid), myByHi(myThid)
289     DO bi = myBxLo(myThid), myBxHi(myThid)
290     DO j=1,sNy
291     DO i=1,sNx
292     x = xC(i,j,bi,bj)
293     y = yC(i,j,bi,bj)
294     C_basal_friction(i,j,bi,bj) =
295     & sqrt(C_basal_fric_const**2*
296     & (1+sin(2*streamice_kx_b_init*PI*x/lenx)*
297     & sin(2*streamice_ky_b_init*PI*y/leny)))
298     ENDDO
299     ENDDO
300     ENDDO
301     ENDDO
302    
303     ELSE IF (STREAMICEbasalTracConfig.EQ.'1DPERIODIC' ) THEN
304    
305     lenx = sNx*nSx*nPx*delX(1)
306     print *, 'lenx', lenx
307     DO bj = myByLo(myThid), myByHi(myThid)
308     DO bi = myBxLo(myThid), myBxHi(myThid)
309     DO j=1,sNy
310     DO i=1,sNx
311     x = xC(i,j,bi,bj)
312     y = yC(i,j,bi,bj)
313     C_basal_friction(i,j,bi,bj) =
314     & sqrt(C_basal_fric_const**2*(1+
315     & sin(2*streamice_kx_b_init*PI*x/lenx)))
316     ENDDO
317     ENDDO
318     ENDDO
319     ENDDO
320    
321     ELSE
322    
323 dgoldberg 1.8 WRITE(msgBuf,'(A)') 'INIT TRAC - NOT IMPLENTED'
324     CALL PRINT_ERROR( msgBuf, myThid)
325     STOP 'ABNORMAL END: S/R STREAMICE_INIT_VAR'
326 dgoldberg 1.4 ENDIF
327    
328     ! finish initialize basal traction
329    
330     CALL STREAMICE_UPD_FFRAC_UNCOUPLED ( myThid )
331 heimbach 1.1
332     _EXCH_XY_RL(H_streamice, myThid )
333     _EXCH_XY_RL(STREAMICE_hmask, myThid )
334     _EXCH_XY_RL(area_shelf_streamice, myThid )
335 dgoldberg 1.4 _EXCH_XY_RL(C_basal_friction, myThid )
336 dgoldberg 1.6 #ifdef USE_ALT_RLOW
337     _EXCH_XY_RL(R_low_si, myThid )
338     #endif
339 dgoldberg 1.4
340 dgoldberg 1.9 !#ifdef STREAMICE_HYBRID_STRESS
341    
342     ! CALL STREAMICE_VISC_BETA (myThid)
343    
344     ! DNG THIS CALL IS TO INITIALISE VISCOSITY
345     ! TO AVOID POSSIBLE ADJOINT INSTABILITIES
346     ! IT IS WRITTEN OVER IN FIRST TIMESTEP
347    
348     #ifdef ALLOW_AUTODIFF
349 dgoldberg 1.4
350 dgoldberg 1.9 CALL STREAMICE_UPD_FFRAC_UNCOUPLED ( myThid )
351     CALL STREAMICE_VELMASK_UPD (myThid)
352     CALL STREAMICE_VEL_SOLVE( myThid )
353 dgoldberg 1.4
354     #endif
355 dgoldberg 1.9
356     !#endif
357 dgoldberg 1.6
358     CALL WRITE_FLD_XY_RL ( "C_basal_fric", "",
359     & C_basal_friction, 0, myThid )
360 heimbach 1.1 CALL WRITE_FLD_XY_RL ( "H_streamice", "init",
361     & H_streamIce, 0, myThid )
362     CALL WRITE_FLD_XY_RL ( "area_shelf_streamice", "init",
363     & area_shelf_streamice, 0, myThid )
364     CALL WRITE_FLD_XY_RL ( "STREAMICE_hmask", "init",
365     & STREAMICE_hmask, 0, myThid )
366 dgoldberg 1.9 #ifdef ALLOW_CTRL
367     CALL ACTIVE_WRITE_GEN_RS( 'maskCtrlst', STREAMICE_ctrl_mask,
368     & 'XY', Nr, 1, .FALSE., 0, mythid, dummyRS )
369     #endif
370     ! call active_write_xyz( 'maskCtrlS', STREAMICE_ctrl_mask, 1, 0,
371     ! & mythid, dummy)
372 dgoldberg 1.2 ! CALL STREAMICE_VELMASK_UPD (myThid)
373     ! CALL STREAMICE_UPD_FFRAC_UNCOUPLED ( myThid )
374     ! CALL STREAMICE_VEL_SOLVE( myThid )
375 heimbach 1.1
376     CALL WRITE_FLD_XY_RL ( "U_init", "",
377 dgoldberg 1.4 & C_basal_friction, 0, myThid )
378 heimbach 1.1 CALL WRITE_FLD_XY_RL ( "V_init", "",
379     & V_streamice, 0, myThid )
380 dgoldberg 1.6 #ifdef USE_ALT_RLOW
381     CALL WRITE_FLD_XY_RL ( "R_low_si", "init",
382     & R_low_si, 0, myThid )
383     #endif
384 heimbach 1.1
385     ! CALL WRITE_FULLARRAY_RL ("H",H_streamice,1,0,0,1,0,myThid)
386     ! CALL WRITE_FULLARRAY_RL ("hmask",STREAMICE_hmask,1,0,0,1,0,myThid)
387     ! CALL WRITE_FULLARRAY_RL ("umask",STREAMICE_umask,1,0,0,1,0,myThid)
388    
389     #endif /* ALLOW_STREAMICE */
390    
391     RETURN
392     END
393    

  ViewVC Help
Powered by ViewVC 1.1.22