dgoldberg/streamice/streamice_init_varia.F

C $Header: /u/gcmpack/MITgcm_contrib/dgoldberg/streamice/streamice_init_varia.F,v 1.15 2013/06/11 17:42:17 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) = streamice_adot_uniform
          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_STREAMICE_2DTRACER
#ifdef STREAMICE_TRACER_AB
         GAD_trac_2d (i,j,bi,bj) = 0. _d 0
!         GAD_trac_2dNm1 (i,j,bi,bj) = 0. _d 0
#endif
#endif
#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)
       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)
       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 

#ifdef ALLOW_STREAMICE_2DTRACER

      IF ( STREAMICETRAC2DINITFILE .NE. ' ' ) THEN
         _BARRIER

        CALL READ_FLD_XY_RL( STREAMICETRAC2dInitFile, ' ',
     &      trac2d, 0, myThid )

       ELSE
        WRITE(msgBuf,'(A)') 'TRAC2dInit - NO FILE SPECIFIED'
        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                    SQUEEZE_RIGHT , 1)
        DO bj = myByLo(myThid), myByHi(myThid)
         DO bi = myBxLo(myThid), myBxHi(myThid)
          DO j=1,sNy
           DO i=1,sNx
            trac2d(i,j,bi,bj) = 0.0
           ENDDO
          ENDDO
         ENDDO
        ENDDO

       ENDIF

      _EXCH_XY_RL (trac2d, myThid)


#endif /*STREAMICE_ALLOW_2DTRACER*/

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