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

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

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

revision 1.3 by dgoldberg, Tue Sep 4 21:11:44 2012 UTC revision 1.12 by dgoldberg, Thu Mar 7 15:23:19 2013 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    #ifdef ALLOW_COST
4    # include "COST_OPTIONS.h"
5    #endif
6  #include "STREAMICE_OPTIONS.h"  #include "STREAMICE_OPTIONS.h"
7    
8  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
# Line 18  C     \================================= Line 20  C     \=================================
20  C     === Global variables ===  C     === Global variables ===
21  #include "SIZE.h"  #include "SIZE.h"
22  #include "GRID.h"  #include "GRID.h"
23    #include "SET_GRID.h"
24  #include "EEPARAMS.h"  #include "EEPARAMS.h"
25  #include "PARAMS.h"  #include "PARAMS.h"
26  #include "STREAMICE.h"  #include "STREAMICE.h"
# Line 34  C     === Local variables === Line 37  C     === Local variables ===
37  C     I,J,bi,bj - Loop counters  C     I,J,bi,bj - Loop counters
38        INTEGER i, j, k, bi, bj, Gi, Gj        INTEGER i, j, k, bi, bj, Gi, Gj
39        INTEGER col_y, col_x        INTEGER col_y, col_x
40        _RL slope_pos, c1        _RL slope_pos, c1, x, y, lenx, leny
41        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
42          _RS     dummyRS
43    
44  CEOP  CEOP
45    
46  C     ZERO OUT FLOATING POINT ARRAYS  C     ZERO OUT FLOATING POINT ARRAYS
47    
48    
49        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
50         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
51          DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
# Line 55  C     ZERO OUT FLOATING POINT ARRAYS Line 61  C     ZERO OUT FLOATING POINT ARRAYS
61            area_shelf_streamice(i,j,bi,bj) = 0. _d 0            area_shelf_streamice(i,j,bi,bj) = 0. _d 0
62            mass_ice_streamice(i,j,bi,bj) = 0. _d 0            mass_ice_streamice(i,j,bi,bj) = 0. _d 0
63            BDOT_streamice(i,j,bi,bj) = 0. _d 0            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            C_basal_friction(i,j,bi,bj) = C_basal_fric_const
66            A_glen(i,j,bi,bj) = A_glen_isothermal            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  #ifdef ALLOW_AUTODIFF_TAMC
69            ru_old_si(i,j,bi,bj) = 0. _d 0            ru_old_si(i,j,bi,bj) = 0. _d 0
70            rv_old_si(i,j,bi,bj) = 0. _d 0            rv_old_si(i,j,bi,bj) = 0. _d 0
# Line 64  C     ZERO OUT FLOATING POINT ARRAYS Line 72  C     ZERO OUT FLOATING POINT ARRAYS
72            zv_old_si(i,j,bi,bj) = 0. _d 0            zv_old_si(i,j,bi,bj) = 0. _d 0
73            h_after_uflux_SI(i,j,bi,bj) = 0. _d 0            h_after_uflux_SI(i,j,bi,bj) = 0. _d 0
74  #endif  #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           ENDDO
88          ENDDO          ENDDO
89    
90    #ifdef ALLOW_COST_TEST
91            cost_func1_streamice (bi,bj) = 0.0
92    #endif
93    
94         ENDDO         ENDDO
95        ENDDO        ENDDO
96    
# Line 105  C     INIT. INTEGER ARRAYS Line 130  C     INIT. INTEGER ARRAYS
130         ENDDO         ENDDO
131        ENDDO        ENDDO
132    
133  !ph      SELECT CASE (TRIM(STREAMICEthickInit))  
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        IF ( STREAMICEthickInit.EQ.'PARAM' ) THEN
154    
# Line 147  C             IF (flow_dir .EQ. 2.0) THE Line 190  C             IF (flow_dir .EQ. 2.0) THE
190                 area_shelf_streamice(i,j,bi,bj) = rA(i,j,bi,bj) *                 area_shelf_streamice(i,j,bi,bj) = rA(i,j,bi,bj) *
191       &          (shelf_edge_pos-xG(i,j,bi,bj)) /       &          (shelf_edge_pos-xG(i,j,bi,bj)) /
192       &          (xG(i+1,j,bi,bj)-xG(i,j,bi,bj))       &          (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                 IF (area_shelf_streamice(i,j,bi,bj).gt. 0. _d 0) THEN
194                  STREAMICE_hmask(i,j,bi,bj) = 2.0                  STREAMICE_hmask(i,j,bi,bj) = 2.0
195                 ELSE                 ELSE
196                  STREAMICE_hmask(i,j,bi,bj) = 0.0                  STREAMICE_hmask(i,j,bi,bj) = 0.0
# Line 167  C             IF (flow_dir .EQ. 2.0) THE Line 210  C             IF (flow_dir .EQ. 2.0) THE
210          ENDDO          ENDDO
211         ENDDO         ENDDO
212    
   
   
213        ELSE IF ( STREAMICEthickInit.EQ.'FILE' ) THEN        ELSE IF ( STREAMICEthickInit.EQ.'FILE' ) THEN
214    
215         IF ( STREAMICEthickFile .NE. ' ' ) THEN         IF ( STREAMICEthickFile .NE. ' ' ) THEN
216          _BARRIER          _BARRIER
217  C The 0 is the "iteration" argument. The ' ' is an empty suffix  C The 0 is the "iteration" argument. The ' ' is an empty suffix
218          CALL READ_FLD_XY_RS( STREAMICEthickFile, ' ', H_streamice,          CALL READ_FLD_XY_RL( STREAMICEthickFile, ' ', H_streamice,
219       &      0, myThid )       &      0, myThid )
220          DO bj = myByLo(myThid), myByHi(myThid)          DO bj = myByLo(myThid), myByHi(myThid)
221           DO bi = myBxLo(myThid), myBxHi(myThid)           DO bi = myBxLo(myThid), myBxHi(myThid)
# Line 182  C The 0 is the "iteration" argument. The Line 223  C The 0 is the "iteration" argument. The
223             DO i=1,sNx             DO i=1,sNx
224              Gi = (myXGlobalLo-1)+(bi-1)*sNx+i              Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
225              Gj = (myYGlobalLo-1)+(bj-1)*sNy+j              Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
226              IF ((Gi.lt.Nx).and.(Gj.lt.Ny)) THEN              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               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)                area_shelf_streamice(i,j,bi,bj) = rA(i,j,bi,bj)
230                STREAMICE_hmask(i,j,bi,bj) = 1.0                STREAMICE_hmask(i,j,bi,bj) = 1.0
# Line 190  C The 0 is the "iteration" argument. The Line 232  C The 0 is the "iteration" argument. The
232                area_shelf_streamice(i,j,bi,bj) = 0. _d 0                area_shelf_streamice(i,j,bi,bj) = 0. _d 0
233                STREAMICE_hmask(i,j,bi,bj) = 0. _d 0                STREAMICE_hmask(i,j,bi,bj) = 0. _d 0
234               ENDIF               ENDIF
235                 Do k=1,Nr
236                 STREAMICE_ctrl_mask(i,j,k,bi,bj) = 1. _d 0
237                 enddo
238              ENDIF              ENDIF
239             ENDDO             ENDDO
240            ENDDO            ENDDO
# Line 208  C The 0 is the "iteration" argument. The Line 253  C The 0 is the "iteration" argument. The
253       &                    SQUEEZE_RIGHT , 1)       &                    SQUEEZE_RIGHT , 1)
254        ENDIF        ENDIF
255    
256        CALL STREAMICE_UPD_FFRAC_UNCOUPLED ( myThid )      #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 basal traction
323    
324          IF ( STREAMICEbasalTracConfig.EQ.'FILE' ) THEN
325    
326           IF ( STREAMICEbasalTracFile .NE. ' ' ) THEN
327            _BARRIER
328    
329            CALL READ_FLD_XY_RL( STREAMICEbasalTracFile, ' ',
330         &      C_basal_friction, 0, myThid )
331    
332           ELSE
333            WRITE(msgBuf,'(A)') 'INIT C_BASAL - FILENAME MISSING'
334            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
335         &                    SQUEEZE_RIGHT , 1)
336           ENDIF
337          
338          ELSE IF (STREAMICEbasalTracConfig.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                C_basal_friction(i,j,bi,bj) = C_basal_fric_const
345               ENDDO
346              ENDDO
347             ENDDO
348            ENDDO
349    
350          ELSE IF (STREAMICEbasalTracConfig.EQ.'2DPERIODIC' ) THEN
351                
352           lenx = sNx*nSx*nPx*delX(1)
353           leny = sNy*nSy*nPy*delY(1)
354           print *, 'lenx', lenx
355           print *, 'leny', leny
356           DO bj = myByLo(myThid), myByHi(myThid)
357            DO bi = myBxLo(myThid), myBxHi(myThid)
358             DO j=1,sNy
359              DO i=1,sNx
360               x = xC(i,j,bi,bj)
361               y = yC(i,j,bi,bj)
362               C_basal_friction(i,j,bi,bj) =
363         &      sqrt(C_basal_fric_const**2*
364         &        (1+sin(2*streamice_kx_b_init*PI*x/lenx)*
365         &           sin(2*streamice_ky_b_init*PI*y/leny)))
366              ENDDO
367             ENDDO
368            ENDDO
369           ENDDO
370    
371          ELSE IF (STREAMICEbasalTracConfig.EQ.'1DPERIODIC' ) THEN
372                
373           lenx = sNx*nSx*nPx*delX(1)
374           print *, 'lenx', lenx
375           DO bj = myByLo(myThid), myByHi(myThid)
376            DO bi = myBxLo(myThid), myBxHi(myThid)
377             DO j=1,sNy
378              DO i=1,sNx
379               x = xC(i,j,bi,bj)
380               y = yC(i,j,bi,bj)
381               C_basal_friction(i,j,bi,bj) =
382         &      sqrt(C_basal_fric_const**2*(1+
383         &        sin(2*streamice_kx_b_init*PI*x/lenx)))
384              ENDDO
385             ENDDO
386            ENDDO
387           ENDDO
388    
389          ELSE
390    
391           WRITE(msgBuf,'(A)') 'INIT TRAC - NOT IMPLENTED'
392           CALL PRINT_ERROR( msgBuf, myThid)
393           STOP 'ABNORMAL END: S/R STREAMICE_INIT_VAR'
394          ENDIF
395    
396    ! finish initialize basal trac
397    
398          CALL STREAMICE_UPD_FFRAC_UNCOUPLED ( myThid )
399    
400        _EXCH_XY_RL(H_streamice, myThid )        _EXCH_XY_RL(H_streamice, myThid )
401        _EXCH_XY_RL(STREAMICE_hmask, myThid )        _EXCH_XY_RL(STREAMICE_hmask, myThid )
402        _EXCH_XY_RL(area_shelf_streamice, myThid )        _EXCH_XY_RL(area_shelf_streamice, myThid )
403          _EXCH_XY_RL(C_basal_friction, myThid )
404          _EXCH_XY_RL(B_glen, myThid )
405    #ifdef USE_ALT_RLOW
406          _EXCH_XY_RL(R_low_si, myThid )
407    #endif
408    
409    !#ifdef STREAMICE_HYBRID_STRESS
410    
411    !      CALL STREAMICE_VISC_BETA (myThid)
412    
413    ! DNG THIS CALL IS TO INITIALISE VISCOSITY
414    !     TO AVOID POSSIBLE ADJOINT INSTABILITIES
415    !     IT IS WRITTEN OVER IN FIRST TIMESTEP
416    
417    #ifdef ALLOW_AUTODIFF
418    
419           CALL STREAMICE_UPD_FFRAC_UNCOUPLED ( myThid )
420           CALL STREAMICE_VELMASK_UPD (myThid)
421           CALL STREAMICE_VEL_SOLVE( myThid )
422    
423    #endif
424    
425    !#endif
426          
427          CALL WRITE_FLD_XY_RL ( "C_basal_fric", "",
428         & C_basal_friction, 0, myThid )
429          CALL WRITE_FLD_XY_RL ( "B_glen_sqrt", "",
430         & B_glen, 0, myThid )
431        CALL WRITE_FLD_XY_RL ( "H_streamice", "init",        CALL WRITE_FLD_XY_RL ( "H_streamice", "init",
432       & H_streamIce, 0, myThid )       & H_streamIce, 0, myThid )
433        CALL WRITE_FLD_XY_RL ( "area_shelf_streamice", "init",        CALL WRITE_FLD_XY_RL ( "area_shelf_streamice", "init",
434       & area_shelf_streamice, 0, myThid )       & area_shelf_streamice, 0, myThid )
435        CALL WRITE_FLD_XY_RL ( "STREAMICE_hmask", "init",        CALL WRITE_FLD_XY_RL ( "STREAMICE_hmask", "init",
436       & STREAMICE_hmask, 0, myThid )       & STREAMICE_hmask, 0, myThid )
437    #ifdef ALLOW_CTRL
438          CALL ACTIVE_WRITE_GEN_RS( 'maskCtrlst', STREAMICE_ctrl_mask,
439         &  'XY', Nr, 1, .FALSE., 0, mythid, dummyRS )
440    #endif
441    !      call active_write_xyz( 'maskCtrlS', STREAMICE_ctrl_mask, 1, 0,
442    !     & mythid, dummy)
443  !       CALL STREAMICE_VELMASK_UPD (myThid)  !       CALL STREAMICE_VELMASK_UPD (myThid)
444  !       CALL STREAMICE_UPD_FFRAC_UNCOUPLED ( myThid )  !       CALL STREAMICE_UPD_FFRAC_UNCOUPLED ( myThid )
445  !       CALL STREAMICE_VEL_SOLVE( myThid )  !       CALL STREAMICE_VEL_SOLVE( myThid )
446    
447        CALL WRITE_FLD_XY_RL ( "U_init", "",        CALL WRITE_FLD_XY_RL ( "U_init", "",
448       &   U_streamice, 0, myThid )       &   C_basal_friction, 0, myThid )
449        CALL WRITE_FLD_XY_RL ( "V_init", "",        CALL WRITE_FLD_XY_RL ( "V_init", "",
450       &   V_streamice, 0, myThid )       &   V_streamice, 0, myThid )
451    #ifdef USE_ALT_RLOW
452          CALL WRITE_FLD_XY_RL ( "R_low_si", "init",
453         & R_low_si, 0, myThid )
454    #endif
455    
456  !       CALL WRITE_FULLARRAY_RL ("H",H_streamice,1,0,0,1,0,myThid)  !       CALL WRITE_FULLARRAY_RL ("H",H_streamice,1,0,0,1,0,myThid)
457  !       CALL WRITE_FULLARRAY_RL ("hmask",STREAMICE_hmask,1,0,0,1,0,myThid)  !       CALL WRITE_FULLARRAY_RL ("hmask",STREAMICE_hmask,1,0,0,1,0,myThid)
458  !       CALL WRITE_FULLARRAY_RL ("umask",STREAMICE_umask,1,0,0,1,0,myThid)  !       CALL WRITE_FULLARRAY_RL ("umask",STREAMICE_umask,1,0,0,1,0,myThid)
459    
460    
461                
462  #endif /* ALLOW_STREAMICE */  #endif /* ALLOW_STREAMICE */
463    
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.12
  ViewVC Help
Powered by ViewVC 1.1.22