/[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.5 - (hide 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 heimbach 1.5 C $Header: /u/gcmpack/MITgcm_contrib/dgoldberg/streamice/streamice_init_varia.F,v 1.4 2012/09/18 17:06:48 dgoldberg Exp $
2 heimbach 1.1 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 dgoldberg 1.4 #include "SET_GRID.h"
22 heimbach 1.1 #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 dgoldberg 1.4 _RL slope_pos, c1, x, y, lenx, leny
39 heimbach 1.1 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 dgoldberg 1.3 BDOT_streamice(i,j,bi,bj) = 0. _d 0
59 dgoldberg 1.4 ! C_basal_friction(i,j,bi,bj) = C_basal_fric_const
60 dgoldberg 1.3 A_glen(i,j,bi,bj) = A_glen_isothermal
61 heimbach 1.1 #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 dgoldberg 1.4 #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 heimbach 1.5 streamice_taubx (i,j,bi,bj) = 0. _d 0
74     streamice_tauby (i,j,bi,bj) = 0. _d 0
75 dgoldberg 1.4 #endif
76 heimbach 1.1 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 dgoldberg 1.4 ! initialize thickness
118 heimbach 1.1
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 heimbach 1.5 IF (area_shelf_streamice(i,j,bi,bj).gt. 0. _d 0) THEN
160 heimbach 1.1 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 dgoldberg 1.4 IF ((Gi.lt.Nx.OR.STREAMICE_EW_periodic).and.
193     & (Gj.lt.Ny.OR.STREAMICE_NS_periodic)) THEN
194 heimbach 1.1 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 dgoldberg 1.4 ! 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 heimbach 1.1
299     _EXCH_XY_RL(H_streamice, myThid )
300     _EXCH_XY_RL(STREAMICE_hmask, myThid )
301     _EXCH_XY_RL(area_shelf_streamice, myThid )
302 dgoldberg 1.4 _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 heimbach 1.1
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 dgoldberg 1.2 ! CALL STREAMICE_VELMASK_UPD (myThid)
319     ! CALL STREAMICE_UPD_FFRAC_UNCOUPLED ( myThid )
320     ! CALL STREAMICE_VEL_SOLVE( myThid )
321 heimbach 1.1
322     CALL WRITE_FLD_XY_RL ( "U_init", "",
323 dgoldberg 1.4 & C_basal_friction, 0, myThid )
324 heimbach 1.1 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