dgoldberg/streamice/streamice_init_varia.F

C $Header: /u/gcmpack/MITgcm_contrib/dgoldberg/streamice/streamice_init_varia.F,v 1.11 2013/01/27 21:05:50 dgoldberg Exp $
C $Name:  $
#ifdef ALLOW_COST
# include "COST_OPTIONS.h"
#endif
#include "STREAMICE_OPTIONS.h"

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

CBOP
      SUBROUTINE STREAMICE_INIT_VARIA( myThid )
C     /============================================================\
C     | SUBROUTINE STREAMICE_INIT_VARIA                             |
C     | o Routine to initialize STREAMICE variables.                |
C     |============================================================|
C     | Initialize STREAMICE parameters and variables.              |
C     \============================================================/
      IMPLICIT NONE

C     === Global variables ===
#include "SIZE.h"
#include "GRID.h"
#include "SET_GRID.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "STREAMICE.h"
#include "STREAMICE_CG.h"
#include "STREAMICE_ADV.h"

C     === Routine arguments ===
C     myThid -  Number of this instance of STREAMICE_INIT_VARIA
      INTEGER myThid
CEndOfInterface

#ifdef ALLOW_STREAMICE
C     === Local variables ===
C     I,J,bi,bj - Loop counters
      INTEGER i, j, k, bi, bj, Gi, Gj
      INTEGER col_y, col_x
      _RL slope_pos, c1, x, y, lenx, leny
      CHARACTER*(MAX_LEN_MBUF) msgBuf
      _RS     dummyRS

CEOP

C     ZERO OUT FLOATING POINT ARRAYS


      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO j=1-Oly,sNy+Oly
         DO i=1-Olx,sNx+Olx
          H_streamIce(i,j,bi,bj) = 0. _d 0
          U_streamice(i,j,bi,bj) = 0. _d 0
          V_streamice(i,j,bi,bj) = 0. _d 0
          visc_streamice(i,j,bi,bj) = 0. _d 0
          tau_beta_eff_streamice(i,j,bi,bj) = 0. _d 0
          float_frac_streamice(i,j,bi,bj) = 0. _d 0
          base_el_streamice(i,j,bi,bj) = 0. _d 0
          surf_el_streamice(i,j,bi,bj) = 0. _d 0
          area_shelf_streamice(i,j,bi,bj) = 0. _d 0
          mass_ice_streamice(i,j,bi,bj) = 0. _d 0
          BDOT_streamice(i,j,bi,bj) = 0. _d 0
          ADOT_streamice(i,j,bi,bj) = 0. _d 0
          C_basal_friction(i,j,bi,bj) = C_basal_fric_const
          B_glen(i,j,bi,bj) = B_glen_isothermal
          H_streamice_prev(i,j,bi,bj) = 0. _d 0
#ifdef ALLOW_AUTODIFF_TAMC
          ru_old_si(i,j,bi,bj) = 0. _d 0
          rv_old_si(i,j,bi,bj) = 0. _d 0
          zu_old_si(i,j,bi,bj) = 0. _d 0
          zv_old_si(i,j,bi,bj) = 0. _d 0
          h_after_uflux_SI(i,j,bi,bj) = 0. _d 0
#endif
#ifdef USE_ALT_RLOW
          R_low_si(i,j,bi,bj) = 0. _d 0
#endif

#ifdef STREAMICE_HYBRID_STRESS
          do k=1,Nr
           visc_streamice_full(i,j,k,bi,bj) =
     &     eps_glen_min**((1-n_glen)/n_glen)
          enddo      
          streamice_taubx (i,j,bi,bj) = 0. _d 0
          streamice_tauby (i,j,bi,bj) = 0. _d 0
#endif
         ENDDO
        ENDDO

#ifdef ALLOW_COST_TEST
        cost_func1_streamice (bi,bj) = 0.0
#endif

       ENDDO
      ENDDO

      DO j = 1-oly, sNy+oly
       DO i = 1-olx, sNx+olx
        DO bj = myByLo(myThid), myByHi(myThid)
         DO bi = myBxLo(myThid), myBxHi(myThid)
cc          DO k=1,4
           DO col_x=-1,1
            DO col_y=-1,1
             streamice_cg_A1(i,j,bi,bj,col_x,col_y)=0.0
             streamice_cg_A2(i,j,bi,bj,col_x,col_y)=0.0
             streamice_cg_A3(i,j,bi,bj,col_x,col_y)=0.0
             streamice_cg_A4(i,j,bi,bj,col_x,col_y)=0.0
            ENDDO
           ENDDO
cc          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO

C     INIT. INTEGER ARRAYS

      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO j=1-Oly,sNy+Oly
         DO i=1-Olx,sNx+Olx
          STREAMICE_hmask(i,j,bi,bj) = -1.0
          STREAMICE_umask(i,j,bi,bj) = 0.0
          STREAMICE_vmask(i,j,bi,bj) = 0.0
          STREAMICE_ufacemask(i,j,bi,bj) = 0.0
          STREAMICE_vfacemask(i,j,bi,bj) = 0.0
          STREAMICE_float_cond(i,j,bi,bj) = 0.0
         ENDDO
        ENDDO
       ENDDO
      ENDDO


#ifdef USE_ALT_RLOW
! init alternate array for topog
      IF ( STREAMICEtopogFile .NE. ' ' ) THEN
        _BARRIER
C The 0 is the "iteration" argument. The ' ' is an empty suffix
       CALL READ_FLD_XY_RS( STREAMICEtopogFile, '', 
     &      R_low_si, 0, myThid )
       
      ELSE
        WRITE(msgBuf,'(A)') 'STREAMICE TOPOG - FILENAME MISSING'
        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                    SQUEEZE_RIGHT , 1)
      ENDIF
#endif

! initialize thickness

#ifndef STREAMICE_GEOM_FILE_SETUP
    
      IF ( STREAMICEthickInit.EQ.'PARAM' ) THEN

      WRITE(msgBuf,'(A)') 'initializing analytic thickness'
      CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                    SQUEEZE_RIGHT , 1)

       slope_pos = shelf_edge_pos - shelf_flat_width
       c1 = 0.0 
       IF (shelf_slope_scale .GT. 0.0) THEN
        c1 = 1.0 / shelf_slope_scale
       ENDIF

       DO bj = myByLo(myThid), myByHi(myThid)
        DO bi = myBxLo(myThid), myBxHi(myThid)
         DO j=1,sNy
          DO i=1,sNx
           Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
           Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
           
           IF ((Gi.lt.Nx).and.(Gj.lt.Ny)) THEN
            
C             IF (flow_dir .EQ. 2.0) THEN
            IF (.TRUE.) THEN
             IF (xC(i-1,j,bi,bj).GE.shelf_edge_pos) THEN
              area_shelf_streamice(i,j,bi,bj) = 0. _d 0
              STREAMICE_hmask(i,j,bi,bj) = 0. _d 0
             ELSE

              IF (xC(i,j,bi,bj).GT.slope_pos) THEN
               H_streamice (i,j,bi,bj) = shelf_min_draft
              ELSE
               H_streamice (i,j,bi,bj) = (shelf_min_draft + 
     &          (shelf_max_draft - shelf_min_draft) * 
     &          min (1.0, (c1*(slope_pos-xC(i,j,bi,bj)))**2))
              ENDIF

              IF (xC(i,j,bi,bj).GT.shelf_edge_pos) THEN
               area_shelf_streamice(i,j,bi,bj) = rA(i,j,bi,bj) *
     &          (shelf_edge_pos-xG(i,j,bi,bj)) /
     &          (xG(i+1,j,bi,bj)-xG(i,j,bi,bj))
               IF (area_shelf_streamice(i,j,bi,bj).gt. 0. _d 0) THEN
                STREAMICE_hmask(i,j,bi,bj) = 2.0
               ELSE
                STREAMICE_hmask(i,j,bi,bj) = 0.0
                H_streamice(i,j,bi,bj) = 0.0
               ENDIF
              ELSE
               area_shelf_streamice(i,j,bi,bj) = rA(i,j,bi,bj) 
               STREAMICE_hmask(i,j,bi,bj) = 1.0
              ENDIF


             ENDIF
            ENDIF
           ENDIF  
          ENDDO
         ENDDO
        ENDDO
       ENDDO

      ELSE IF ( STREAMICEthickInit.EQ.'FILE' ) THEN

       IF ( STREAMICEthickFile .NE. ' ' ) THEN
        _BARRIER
C The 0 is the "iteration" argument. The ' ' is an empty suffix
        CALL READ_FLD_XY_RL( STREAMICEthickFile, ' ', H_streamice,
     &      0, myThid )
        DO bj = myByLo(myThid), myByHi(myThid)
         DO bi = myBxLo(myThid), myBxHi(myThid)
          DO j=1,sNy
           DO i=1,sNx
            Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
            Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
            IF ((Gi.lt.Nx.OR.STREAMICE_EW_periodic).and.
     &          (Gj.lt.Ny.OR.STREAMICE_NS_periodic)) THEN
             IF (H_streamice(i,j,bi,bj).GT.0. _d 0) THEN
              area_shelf_streamice(i,j,bi,bj) = rA(i,j,bi,bj)
              STREAMICE_hmask(i,j,bi,bj) = 1.0
             ELSE
              area_shelf_streamice(i,j,bi,bj) = 0. _d 0
              STREAMICE_hmask(i,j,bi,bj) = 0. _d 0
             ENDIF
             Do k=1,Nr
             STREAMICE_ctrl_mask(i,j,k,bi,bj) = 1. _d 0
             enddo
            ENDIF
           ENDDO
          ENDDO
         ENDDO
        ENDDO
       ELSE
        WRITE(msgBuf,'(A)') 'INIT THICKNESS - FILENAME MISSING'
        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                    SQUEEZE_RIGHT , 1)
       ENDIF
      
      ELSE

       WRITE(msgBuf,'(A)') 'INIT THICKNESS - NOT IMPLENTED'
       CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                    SQUEEZE_RIGHT , 1)
      ENDIF

#else
! STREAMICE_GEOM_FILE_SETUP - init thickness and hmask MUST come from file

      IF ( STREAMICEthickFile .NE. ' ' ) THEN
        _BARRIER
C The 0 is the "iteration" argument. The ' ' is an empty suffix
      CALL READ_FLD_XY_RL( STREAMICEthickFile, ' ', H_streamice,
     &      0, myThid )
      ELSE
       WRITE(msgBuf,'(A)') 'INIT THICKNESS - FILENAME MISSING'
       CALL PRINT_ERROR( msgBuf, myThid)
      ENDIF

      IF ( STREAMICEhMaskFile .NE. ' ' ) THEN
        _BARRIER
C The 0 is the "iteration" argument. The ' ' is an empty suffix
      CALL READ_FLD_XY_RS( STREAMICEhMaskFile, ' ', STREAMICE_hmask,
     &      0, myThid )
      ELSE
       WRITE(msgBuf,'(A)') 'INIT HMASK - FILENAME MISSING'
       CALL PRINT_ERROR( msgBuf, myThid)
      ENDIF

#endif 
! STREAMICE_GEOM_FILE_SETUP


! finish initialize thickness

! initialize glen constant

      IF ( STREAMICEGlenConstConfig.EQ.'FILE' ) THEN

       IF ( STREAMICEGlenConstFile .NE. ' ' ) THEN
        _BARRIER

        CALL READ_FLD_XY_RL( STREAMICEGlenConstFile, ' ', 
     &      B_glen, 0, myThid )

       ELSE
        WRITE(msgBuf,'(A)') 'INIT GLEN - FILENAME MISSING'
        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                    SQUEEZE_RIGHT , 1)
       ENDIF
      
      ELSE IF (STREAMICEGlenConstConfig.EQ.'UNIFORM' ) THEN

        DO bj = myByLo(myThid), myByHi(myThid)
         DO bi = myBxLo(myThid), myBxHi(myThid)
          DO j=1,sNy
           DO i=1,sNx
            B_glen(i,j,bi,bj) = B_glen_isothermal
           ENDDO
          ENDDO
         ENDDO
        ENDDO

      ELSE

       WRITE(msgBuf,'(A)') 'INIT GLEN CONSTANT - NOT IMPLENTED'
       CALL PRINT_ERROR( msgBuf, myThid)
       STOP 'ABNORMAL END: S/R STREAMICE_INIT_VAR'
      ENDIF

! finish initialize glen constant

! initialize basal traction

      IF ( STREAMICEbasalTracConfig.EQ.'FILE' ) THEN

       IF ( STREAMICEbasalTracFile .NE. ' ' ) THEN
        _BARRIER

        CALL READ_FLD_XY_RL( STREAMICEbasalTracFile, ' ', 
     &      C_basal_friction, 0, myThid )

       ELSE
        WRITE(msgBuf,'(A)') 'INIT C_BASAL - FILENAME MISSING'
        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                    SQUEEZE_RIGHT , 1)
       ENDIF
      
      ELSE IF (STREAMICEbasalTracConfig.EQ.'UNIFORM' ) THEN

        DO bj = myByLo(myThid), myByHi(myThid)
         DO bi = myBxLo(myThid), myBxHi(myThid)
          DO j=1,sNy
           DO i=1,sNx
            C_basal_friction(i,j,bi,bj) = C_basal_fric_const
           ENDDO
          ENDDO
         ENDDO
        ENDDO

      ELSE IF (STREAMICEbasalTracConfig.EQ.'2DPERIODIC' ) THEN
            
       lenx = sNx*nSx*nPx*delX(1)
       leny = sNy*nSy*nPy*delY(1)
       print *, 'lenx', lenx
       print *, 'leny', leny
       DO bj = myByLo(myThid), myByHi(myThid)
        DO bi = myBxLo(myThid), myBxHi(myThid)
         DO j=1,sNy
          DO i=1,sNx
           x = xC(i,j,bi,bj)
           y = yC(i,j,bi,bj)
           C_basal_friction(i,j,bi,bj) = 
     &      sqrt(C_basal_fric_const**2*
     &        (1+sin(2*streamice_kx_b_init*PI*x/lenx)*
     &           sin(2*streamice_ky_b_init*PI*y/leny)))
          ENDDO
         ENDDO
        ENDDO
       ENDDO

      ELSE IF (STREAMICEbasalTracConfig.EQ.'1DPERIODIC' ) THEN
            
       lenx = sNx*nSx*nPx*delX(1)
       print *, 'lenx', lenx
       DO bj = myByLo(myThid), myByHi(myThid)
        DO bi = myBxLo(myThid), myBxHi(myThid)
         DO j=1,sNy
          DO i=1,sNx
           x = xC(i,j,bi,bj)
           y = yC(i,j,bi,bj)
           C_basal_friction(i,j,bi,bj) = 
     &      sqrt(C_basal_fric_const**2*(1+
     &        sin(2*streamice_kx_b_init*PI*x/lenx)))
          ENDDO
         ENDDO
        ENDDO
       ENDDO

      ELSE

       WRITE(msgBuf,'(A)') 'INIT TRAC - NOT IMPLENTED'
       CALL PRINT_ERROR( msgBuf, myThid)
       STOP 'ABNORMAL END: S/R STREAMICE_INIT_VAR'
      ENDIF

! finish initialize basal trac 

      CALL STREAMICE_UPD_FFRAC_UNCOUPLED ( myThid )

      _EXCH_XY_RL(H_streamice, myThid )
      _EXCH_XY_RL(STREAMICE_hmask, myThid )
      _EXCH_XY_RL(area_shelf_streamice, myThid )
      _EXCH_XY_RL(C_basal_friction, myThid )
      _EXCH_XY_RL(B_glen, myThid )
#ifdef USE_ALT_RLOW
      _EXCH_XY_RL(R_low_si, myThid )
#endif

!#ifdef STREAMICE_HYBRID_STRESS

!      CALL STREAMICE_VISC_BETA (myThid)

! DNG THIS CALL IS TO INITIALISE VISCOSITY
!     TO AVOID POSSIBLE ADJOINT INSTABILITIES
!     IT IS WRITTEN OVER IN FIRST TIMESTEP

#ifdef ALLOW_AUTODIFF

       CALL STREAMICE_UPD_FFRAC_UNCOUPLED ( myThid )
       CALL STREAMICE_VELMASK_UPD (myThid)
       CALL STREAMICE_VEL_SOLVE( myThid )

#endif

!#endif
      
      CALL WRITE_FLD_XY_RL ( "C_basal_fric", "", 
     & C_basal_friction, 0, myThid )
      CALL WRITE_FLD_XY_RL ( "B_glen_sqrt", "",
     & B_glen, 0, myThid )
      CALL WRITE_FLD_XY_RL ( "H_streamice", "init", 
     & H_streamIce, 0, myThid )
      CALL WRITE_FLD_XY_RL ( "area_shelf_streamice", "init", 
     & area_shelf_streamice, 0, myThid )
      CALL WRITE_FLD_XY_RL ( "STREAMICE_hmask", "init", 
     & STREAMICE_hmask, 0, myThid )
#ifdef ALLOW_CTRL
      CALL ACTIVE_WRITE_GEN_RS( 'maskCtrlst', STREAMICE_ctrl_mask,
     &  'XY', Nr, 1, .FALSE., 0, mythid, dummyRS )
#endif
!      call active_write_xyz( 'maskCtrlS', STREAMICE_ctrl_mask, 1, 0,
!     & mythid, dummy)
!       CALL STREAMICE_VELMASK_UPD (myThid)
!       CALL STREAMICE_UPD_FFRAC_UNCOUPLED ( myThid )
!       CALL STREAMICE_VEL_SOLVE( myThid )

      CALL WRITE_FLD_XY_RL ( "U_init", "", 
     &   C_basal_friction, 0, myThid )
      CALL WRITE_FLD_XY_RL ( "V_init", "", 
     &   V_streamice, 0, myThid )
#ifdef USE_ALT_RLOW
      CALL WRITE_FLD_XY_RL ( "R_low_si", "init", 
     & R_low_si, 0, myThid )
#endif

!       CALL WRITE_FULLARRAY_RL ("H",H_streamice,1,0,0,1,0,myThid)
!       CALL WRITE_FULLARRAY_RL ("hmask",STREAMICE_hmask,1,0,0,1,0,myThid)
!       CALL WRITE_FULLARRAY_RL ("umask",STREAMICE_umask,1,0,0,1,0,myThid)


#endif /* ALLOW_STREAMICE */

      RETURN
      END

1	dgoldberg	1.12	C $Header: /u/gcmpack/MITgcm_contrib/dgoldberg/streamice/streamice_init_varia.F,v 1.11 2013/01/27 21:05:50 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	dgoldberg	1.12
49	heimbach	1.1	DO bj = myByLo(myThid), myByHi(myThid)
50			DO bi = myBxLo(myThid), myBxHi(myThid)
51			DO j=1-Oly,sNy+Oly
52			DO i=1-Olx,sNx+Olx
53			H_streamIce(i,j,bi,bj) = 0. _d 0
54			U_streamice(i,j,bi,bj) = 0. _d 0
55			V_streamice(i,j,bi,bj) = 0. _d 0
56			visc_streamice(i,j,bi,bj) = 0. _d 0
57			tau_beta_eff_streamice(i,j,bi,bj) = 0. _d 0
58			float_frac_streamice(i,j,bi,bj) = 0. _d 0
59			base_el_streamice(i,j,bi,bj) = 0. _d 0
60			surf_el_streamice(i,j,bi,bj) = 0. _d 0
61			area_shelf_streamice(i,j,bi,bj) = 0. _d 0
62			mass_ice_streamice(i,j,bi,bj) = 0. _d 0
63	dgoldberg	1.3	BDOT_streamice(i,j,bi,bj) = 0. _d 0
64	dgoldberg	1.10	ADOT_streamice(i,j,bi,bj) = 0. _d 0
65	dgoldberg	1.12	C_basal_friction(i,j,bi,bj) = C_basal_fric_const
66			B_glen(i,j,bi,bj) = B_glen_isothermal
67	dgoldberg	1.7	H_streamice_prev(i,j,bi,bj) = 0. _d 0
68	heimbach	1.1	#ifdef ALLOW_AUTODIFF_TAMC
69			ru_old_si(i,j,bi,bj) = 0. _d 0
70			rv_old_si(i,j,bi,bj) = 0. _d 0
71			zu_old_si(i,j,bi,bj) = 0. _d 0
72			zv_old_si(i,j,bi,bj) = 0. _d 0
73			h_after_uflux_SI(i,j,bi,bj) = 0. _d 0
74			#endif
75	dgoldberg	1.6	#ifdef USE_ALT_RLOW
76			R_low_si(i,j,bi,bj) = 0. _d 0
77			#endif
78
79	dgoldberg	1.4	#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	heimbach	1.5	streamice_taubx (i,j,bi,bj) = 0. _d 0
85			streamice_tauby (i,j,bi,bj) = 0. _d 0
86	dgoldberg	1.4	#endif
87	heimbach	1.1	ENDDO
88			ENDDO
89	dgoldberg	1.6
90			#ifdef ALLOW_COST_TEST
91			cost_func1_streamice (bi,bj) = 0.0
92			#endif
93
94	heimbach	1.1	ENDDO
95			ENDDO
96
97			DO j = 1-oly, sNy+oly
98			DO i = 1-olx, sNx+olx
99			DO bj = myByLo(myThid), myByHi(myThid)
100			DO bi = myBxLo(myThid), myBxHi(myThid)
101			cc DO k=1,4
102			DO col_x=-1,1
103			DO col_y=-1,1
104			streamice_cg_A1(i,j,bi,bj,col_x,col_y)=0.0
105			streamice_cg_A2(i,j,bi,bj,col_x,col_y)=0.0
106			streamice_cg_A3(i,j,bi,bj,col_x,col_y)=0.0
107			streamice_cg_A4(i,j,bi,bj,col_x,col_y)=0.0
108			ENDDO
109			ENDDO
110			cc ENDDO
111			ENDDO
112			ENDDO
113			ENDDO
114			ENDDO
115
116			C INIT. INTEGER ARRAYS
117
118			DO bj = myByLo(myThid), myByHi(myThid)
119			DO bi = myBxLo(myThid), myBxHi(myThid)
120			DO j=1-Oly,sNy+Oly
121			DO i=1-Olx,sNx+Olx
122			STREAMICE_hmask(i,j,bi,bj) = -1.0
123			STREAMICE_umask(i,j,bi,bj) = 0.0
124			STREAMICE_vmask(i,j,bi,bj) = 0.0
125			STREAMICE_ufacemask(i,j,bi,bj) = 0.0
126			STREAMICE_vfacemask(i,j,bi,bj) = 0.0
127			STREAMICE_float_cond(i,j,bi,bj) = 0.0
128			ENDDO
129			ENDDO
130			ENDDO
131			ENDDO
132
133	dgoldberg	1.6
134			#ifdef USE_ALT_RLOW
135			! init alternate array for topog
136	dgoldberg	1.12	IF ( STREAMICEtopogFile .NE. ' ' ) THEN
137	dgoldberg	1.6	_BARRIER
138			C The 0 is the "iteration" argument. The ' ' is an empty suffix
139	dgoldberg	1.12	CALL READ_FLD_XY_RS( STREAMICEtopogFile, '',
140	dgoldberg	1.6	& 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	dgoldberg	1.4	! initialize thickness
150	dgoldberg	1.10
151			#ifndef STREAMICE_GEOM_FILE_SETUP
152	heimbach	1.1
153			IF ( STREAMICEthickInit.EQ.'PARAM' ) THEN
154
155			WRITE(msgBuf,'(A)') 'initializing analytic thickness'
156			CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
157			& SQUEEZE_RIGHT , 1)
158
159			slope_pos = shelf_edge_pos - shelf_flat_width
160			c1 = 0.0
161			IF (shelf_slope_scale .GT. 0.0) THEN
162			c1 = 1.0 / shelf_slope_scale
163			ENDIF
164
165			DO bj = myByLo(myThid), myByHi(myThid)
166			DO bi = myBxLo(myThid), myBxHi(myThid)
167			DO j=1,sNy
168			DO i=1,sNx
169			Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
170			Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
171
172			IF ((Gi.lt.Nx).and.(Gj.lt.Ny)) THEN
173
174			C IF (flow_dir .EQ. 2.0) THEN
175			IF (.TRUE.) THEN
176			IF (xC(i-1,j,bi,bj).GE.shelf_edge_pos) THEN
177			area_shelf_streamice(i,j,bi,bj) = 0. _d 0
178			STREAMICE_hmask(i,j,bi,bj) = 0. _d 0
179			ELSE
180
181			IF (xC(i,j,bi,bj).GT.slope_pos) THEN
182			H_streamice (i,j,bi,bj) = shelf_min_draft
183			ELSE
184			H_streamice (i,j,bi,bj) = (shelf_min_draft +
185			& (shelf_max_draft - shelf_min_draft) *
186			& min (1.0, (c1(slope_pos-xC(i,j,bi,bj)))*2))
187			ENDIF
188
189			IF (xC(i,j,bi,bj).GT.shelf_edge_pos) THEN
190			area_shelf_streamice(i,j,bi,bj) = rA(i,j,bi,bj) *
191			& (shelf_edge_pos-xG(i,j,bi,bj)) /
192			& (xG(i+1,j,bi,bj)-xG(i,j,bi,bj))
193	heimbach	1.5	IF (area_shelf_streamice(i,j,bi,bj).gt. 0. _d 0) THEN
194	heimbach	1.1	STREAMICE_hmask(i,j,bi,bj) = 2.0
195			ELSE
196			STREAMICE_hmask(i,j,bi,bj) = 0.0
197			H_streamice(i,j,bi,bj) = 0.0
198			ENDIF
199			ELSE
200			area_shelf_streamice(i,j,bi,bj) = rA(i,j,bi,bj)
201			STREAMICE_hmask(i,j,bi,bj) = 1.0
202			ENDIF
203
204
205			ENDIF
206			ENDIF
207			ENDIF
208			ENDDO
209			ENDDO
210			ENDDO
211			ENDDO
212
213			ELSE IF ( STREAMICEthickInit.EQ.'FILE' ) THEN
214
215			IF ( STREAMICEthickFile .NE. ' ' ) THEN
216			_BARRIER
217			C The 0 is the "iteration" argument. The ' ' is an empty suffix
218	dgoldberg	1.10	CALL READ_FLD_XY_RL( STREAMICEthickFile, ' ', H_streamice,
219	heimbach	1.1	& 0, myThid )
220			DO bj = myByLo(myThid), myByHi(myThid)
221			DO bi = myBxLo(myThid), myBxHi(myThid)
222			DO j=1,sNy
223			DO i=1,sNx
224			Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
225			Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
226	dgoldberg	1.4	IF ((Gi.lt.Nx.OR.STREAMICE_EW_periodic).and.
227			& (Gj.lt.Ny.OR.STREAMICE_NS_periodic)) THEN
228	heimbach	1.1	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)
230			STREAMICE_hmask(i,j,bi,bj) = 1.0
231			ELSE
232			area_shelf_streamice(i,j,bi,bj) = 0. _d 0
233			STREAMICE_hmask(i,j,bi,bj) = 0. _d 0
234			ENDIF
235	dgoldberg	1.9	Do k=1,Nr
236	dgoldberg	1.12	STREAMICE_ctrl_mask(i,j,k,bi,bj) = 1. _d 0
237	dgoldberg	1.9	enddo
238	heimbach	1.1	ENDIF
239			ENDDO
240			ENDDO
241			ENDDO
242			ENDDO
243			ELSE
244			WRITE(msgBuf,'(A)') 'INIT THICKNESS - FILENAME MISSING'
245			CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
246			& SQUEEZE_RIGHT , 1)
247			ENDIF
248
249			ELSE
250
251			WRITE(msgBuf,'(A)') 'INIT THICKNESS - NOT IMPLENTED'
252			CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
253			& SQUEEZE_RIGHT , 1)
254			ENDIF
255
256	dgoldberg	1.10	#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	dgoldberg	1.4	! finish initialize thickness
284
285	dgoldberg	1.11	! initialize glen constant
286	dgoldberg	1.4
287	dgoldberg	1.11	IF ( STREAMICEGlenConstConfig.EQ.'FILE' ) THEN
288	dgoldberg	1.4
289	dgoldberg	1.11	IF ( STREAMICEGlenConstFile .NE. ' ' ) THEN
290	dgoldberg	1.4	_BARRIER
291	dgoldberg	1.11
292			CALL READ_FLD_XY_RL( STREAMICEGlenConstFile, ' ',
293	dgoldberg	1.12	& B_glen, 0, myThid )
294	dgoldberg	1.4
295			ELSE
296	dgoldberg	1.11	WRITE(msgBuf,'(A)') 'INIT GLEN - FILENAME MISSING'
297	dgoldberg	1.4	CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
298			& SQUEEZE_RIGHT , 1)
299			ENDIF
300
301	dgoldberg	1.11	ELSE IF (STREAMICEGlenConstConfig.EQ.'UNIFORM' ) THEN
302	dgoldberg	1.4
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	dgoldberg	1.12	B_glen(i,j,bi,bj) = B_glen_isothermal
308	dgoldberg	1.10	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	dgoldberg	1.11	! finish initialize glen constant
321
322			! initialize basal traction
323	dgoldberg	1.10
324	dgoldberg	1.11	IF ( STREAMICEbasalTracConfig.EQ.'FILE' ) THEN
325	dgoldberg	1.10
326	dgoldberg	1.11	IF ( STREAMICEbasalTracFile .NE. ' ' ) THEN
327	dgoldberg	1.10	_BARRIER
328	dgoldberg	1.11
329			CALL READ_FLD_XY_RL( STREAMICEbasalTracFile, ' ',
330			& C_basal_friction, 0, myThid )
331	dgoldberg	1.10
332			ELSE
333	dgoldberg	1.11	WRITE(msgBuf,'(A)') 'INIT C_BASAL - FILENAME MISSING'
334	dgoldberg	1.10	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	dgoldberg	1.4	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 = sNxnSxnPx*delX(1)
353			leny = sNynSynPy*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(2streamice_kx_b_initPIx/lenx)
365			& sin(2streamice_ky_b_initPI*y/leny)))
366			ENDDO
367			ENDDO
368			ENDDO
369			ENDDO
370
371			ELSE IF (STREAMICEbasalTracConfig.EQ.'1DPERIODIC' ) THEN
372
373			lenx = sNxnSxnPx*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(2streamice_kx_b_initPI*x/lenx)))
384			ENDDO
385			ENDDO
386			ENDDO
387			ENDDO
388
389			ELSE
390
391	dgoldberg	1.8	WRITE(msgBuf,'(A)') 'INIT TRAC - NOT IMPLENTED'
392			CALL PRINT_ERROR( msgBuf, myThid)
393			STOP 'ABNORMAL END: S/R STREAMICE_INIT_VAR'
394	dgoldberg	1.4	ENDIF
395
396	dgoldberg	1.11	! finish initialize basal trac
397	dgoldberg	1.4
398			CALL STREAMICE_UPD_FFRAC_UNCOUPLED ( myThid )
399	heimbach	1.1
400			_EXCH_XY_RL(H_streamice, myThid )
401			_EXCH_XY_RL(STREAMICE_hmask, myThid )
402			_EXCH_XY_RL(area_shelf_streamice, myThid )
403	dgoldberg	1.4	_EXCH_XY_RL(C_basal_friction, myThid )
404	dgoldberg	1.12	_EXCH_XY_RL(B_glen, myThid )
405	dgoldberg	1.6	#ifdef USE_ALT_RLOW
406			_EXCH_XY_RL(R_low_si, myThid )
407			#endif
408	dgoldberg	1.4
409	dgoldberg	1.9	!#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	dgoldberg	1.4
419	dgoldberg	1.9	CALL STREAMICE_UPD_FFRAC_UNCOUPLED ( myThid )
420			CALL STREAMICE_VELMASK_UPD (myThid)
421			CALL STREAMICE_VEL_SOLVE( myThid )
422	dgoldberg	1.4
423			#endif
424	dgoldberg	1.9
425			!#endif
426	dgoldberg	1.6
427			CALL WRITE_FLD_XY_RL ( "C_basal_fric", "",
428			& C_basal_friction, 0, myThid )
429	dgoldberg	1.12	CALL WRITE_FLD_XY_RL ( "B_glen_sqrt", "",
430			& B_glen, 0, myThid )
431	heimbach	1.1	CALL WRITE_FLD_XY_RL ( "H_streamice", "init",
432			& H_streamIce, 0, myThid )
433			CALL WRITE_FLD_XY_RL ( "area_shelf_streamice", "init",
434			& area_shelf_streamice, 0, myThid )
435			CALL WRITE_FLD_XY_RL ( "STREAMICE_hmask", "init",
436			& STREAMICE_hmask, 0, myThid )
437	dgoldberg	1.9	#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	dgoldberg	1.2	! CALL STREAMICE_VELMASK_UPD (myThid)
444			! CALL STREAMICE_UPD_FFRAC_UNCOUPLED ( myThid )
445			! CALL STREAMICE_VEL_SOLVE( myThid )
446	heimbach	1.1
447			CALL WRITE_FLD_XY_RL ( "U_init", "",
448	dgoldberg	1.4	& C_basal_friction, 0, myThid )
449	heimbach	1.1	CALL WRITE_FLD_XY_RL ( "V_init", "",
450			& V_streamice, 0, myThid )
451	dgoldberg	1.6	#ifdef USE_ALT_RLOW
452			CALL WRITE_FLD_XY_RL ( "R_low_si", "init",
453			& R_low_si, 0, myThid )
454			#endif
455	heimbach	1.1
456			! 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)
458			! CALL WRITE_FULLARRAY_RL ("umask",STREAMICE_umask,1,0,0,1,0,myThid)
459	dgoldberg	1.12
460
461	heimbach	1.1
462			#endif /* ALLOW_STREAMICE */
463
464			RETURN
465			END
466