/[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.11 - (hide annotations) (download)
Sun Jan 27 21:05:50 2013 UTC (12 years, 5 months ago) by dgoldberg
Branch: MAIN
Changes since 1.10: +19 -17 lines
allow separate mask maskCtrls (as opposed to maskCtrlC) for ctrl_pack and unpack

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

  ViewVC Help
Powered by ViewVC 1.1.22