pkg/seaice/seaice_init_varia.F

C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_init_varia.F,v 1.55 2011/04/20 01:34:27 jmc Exp $
C $Name:  $

#include "SEAICE_OPTIONS.h"

CStartOfInterface
      SUBROUTINE SEAICE_INIT_VARIA( myThid )
C     *==========================================================*
C     | SUBROUTINE SEAICE_INIT_VARIA                             |
C     | o Initialization of sea ice model.                       |
C     *==========================================================*
C     *==========================================================*
      IMPLICIT NONE

C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "DYNVARS.h"
#include "FFIELDS.h"
#include "SEAICE_SIZE.h"
#include "SEAICE_PARAMS.h"
#include "SEAICE.h"
#include "SEAICE_TRACER.h"
#include "SEAICE_TAVE.h"
#ifdef ALLOW_EXCH2
# include "W2_EXCH2_SIZE.h"
# include "W2_EXCH2_TOPOLOGY.h"
#endif
#ifdef ALLOW_OBCS
# include "OBCS_OPTIONS.h"
# include "OBCS.h"
#endif

C     === Routine arguments ===
C     myThid - Thread no. that called this routine.
      INTEGER myThid
CEndOfInterface

C     === Local variables ===
C     i,j,k,bi,bj - Loop counters

      INTEGER i, j, bi, bj
      _RL PSTAR
      INTEGER myIter, kSurface
#ifdef SEAICE_CGRID
      _RS  mask_uice
#endif
#ifdef SEAICE_MULTICATEGORY
      INTEGER k
#endif
#ifdef SEAICE_AGE
      INTEGER iTracer
#endif
#ifdef ALLOW_OBCS
      INTEGER I_obc, J_obc
#endif /* ALLOW_OBCS */
#ifdef ALLOW_EXCH2
# ifndef SEAICE_CGRID
      INTEGER myTile
# endif
#endif

      IF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
       kSurface = Nr
      ELSE
       kSurface = 1
      ENDIF

C--   Initialise all variables in common blocks:
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          HEFF(i,j,bi,bj)=0. _d 0
          AREA(i,j,bi,bj)=0. _d 0
          UICE(i,j,bi,bj)=0. _d 0
          VICE(i,j,bi,bj)=0. _d 0
C
          uIceNm1(i,j,bi,bj)=0. _d 0
          vIceNm1(i,j,bi,bj)=0. _d 0
          areaNm1(i,j,bi,bj)=0. _d 0
          hEffNm1(i,j,bi,bj)=0. _d 0
         ENDDO
        ENDDO
#ifdef SEAICE_MULTICATEGORY
        DO k=1,MULTDIM
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           TICES(i,j,k,bi,bj)=0. _d 0
          ENDDO
         ENDDO
        ENDDO
#endif
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          ETA (i,j,bi,bj)   = 0. _d 0
          ZETA(i,j,bi,bj)   = 0. _d 0
          DRAGS(i,j,bi,bj)  = 0. _d 0
          DRAGA(i,j,bi,bj)  = 0. _d 0
          FORCEX(i,j,bi,bj) = 0. _d 0
          FORCEY(i,j,bi,bj) = 0. _d 0
          UICEC(i,j,bi,bj)  = 0. _d 0
          VICEC(i,j,bi,bj)  = 0. _d 0
#ifdef SEAICE_CGRID
          seaiceMassC(i,j,bi,bj)=0. _d 0
          seaiceMassU(i,j,bi,bj)=0. _d 0
          seaiceMassV(i,j,bi,bj)=0. _d 0
          stressDivergenceX(i,j,bi,bj) = 0. _d 0
          stressDivergenceY(i,j,bi,bj) = 0. _d 0
# ifdef SEAICE_ALLOW_EVP
          seaice_sigma1 (i,j,bi,bj) = 0. _d 0
          seaice_sigma2 (i,j,bi,bj) = 0. _d 0
          seaice_sigma12(i,j,bi,bj) = 0. _d 0
# endif /* SEAICE_ALLOW_EVP */
#else /* SEAICE_CGRID */
          AMASS(i,j,bi,bj)  = 0. _d 0
          DAIRN(i,j,bi,bj)  = 0. _d 0
          WINDX(i,j,bi,bj)  = 0. _d 0
          WINDY(i,j,bi,bj)  = 0. _d 0
          GWATX(i,j,bi,bj)  = 0. _d 0
          GWATY(i,j,bi,bj)  = 0. _d 0
#endif /* SEAICE_CGRID */
          DWATN(i,j,bi,bj)  = 0. _d 0
          PRESS0(i,j,bi,bj) = 0. _d 0
          FORCEX0(i,j,bi,bj)= 0. _d 0
          FORCEY0(i,j,bi,bj)= 0. _d 0
          ZMAX(i,j,bi,bj)   = 0. _d 0
          ZMIN(i,j,bi,bj)   = 0. _d 0
          HSNOW(i,j,bi,bj)  = 0. _d 0
#ifdef SEAICE_VARIABLE_SALINITY
          HSALT(i,j,bi,bj)  = 0. _d 0
#endif
#ifdef SEAICE_AGE
          DO iTracer = 1, SEAICE_num
           IceAgeTr(i,j,bi,bj,iTracer) = 0. _d 0
          ENDDO
#endif
          YNEG (i,j,bi,bj)  = 0. _d 0
          RIVER(i,j,bi,bj)  = 0. _d 0
          TMIX(i,j,bi,bj)   = 0. _d 0
          TICE(i,j,bi,bj)   = 0. _d 0
          TAUX(i,j,bi,bj)   = 0. _d 0
          TAUY(i,j,bi,bj)   = 0. _d 0
#ifdef ALLOW_SEAICE_COST_EXPORT
          uHeffExportCell(i,j,bi,bj) = 0. _d 0
          vHeffExportCell(i,j,bi,bj) = 0. _d 0
#endif
          saltWtrIce(i,j,bi,bj) = 0. _d 0
          frWtrIce(i,j,bi,bj)   = 0. _d 0
#ifdef ALLOW_MEAN_SFLUX_COST_CONTRIBUTION
          frWtrAtm(i,j,bi,bj)   = 0. _d 0
#endif
         ENDDO
        ENDDO
       ENDDO
      ENDDO

#ifdef ALLOW_TIMEAVE
C     Initialize averages to zero
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        CALL TIMEAVE_RESET( FUtave   , 1, bi, bj, myThid )
        CALL TIMEAVE_RESET( FVtave   , 1, bi, bj, myThid )
        CALL TIMEAVE_RESET( EmPmRtave, 1, bi, bj, myThid )
        CALL TIMEAVE_RESET( QNETtave , 1, bi, bj, myThid )
        CALL TIMEAVE_RESET( QSWtave  , 1, bi, bj, myThid )
        CALL TIMEAVE_RESET( UICEtave , 1, bi, bj, myThid )
        CALL TIMEAVE_RESET( VICEtave , 1, bi, bj, myThid )
        CALL TIMEAVE_RESET( HEFFtave , 1, bi, bj, myThid )
        CALL TIMEAVE_RESET( AREAtave , 1, bi, bj, myThid )
        SEAICE_timeAve(bi,bj) = ZERO
       ENDDO
      ENDDO
#endif /* ALLOW_TIMEAVE */

C--   Initialize (variable) grid info. As long as we allow masking of
C--   velocities outside of ice covered areas (in seaice_dynsolver)
C--   we need to re-initialize seaiceMaskU/V here for TAF/TAMC
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO j=1-OLy+1,sNy+OLy
         DO i=1-OLx+1,sNx+OLx
#ifdef SEAICE_CGRID
          seaiceMaskU(i,j,bi,bj)=   0.0 _d 0
          seaiceMaskV(i,j,bi,bj)=   0.0 _d 0
          mask_uice=HEFFM(i,j,bi,bj)+HEFFM(i-1,j  ,bi,bj)
          IF(mask_uice.GT.1.5 _d 0) seaiceMaskU(i,j,bi,bj)=1.0 _d 0
          mask_uice=HEFFM(i,j,bi,bj)+HEFFM(i  ,j-1,bi,bj)
          IF(mask_uice.GT.1.5 _d 0) seaiceMaskV(i,j,bi,bj)=1.0 _d 0
#else
          UVM(i,j,bi,bj)=   0.0 _d 0
#endif /* SEAICE_CGRID */
         ENDDO
        ENDDO
       ENDDO
      ENDDO

      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
#ifdef ALLOW_OBCS
        IF (useOBCS) THEN
C--   If OBCS is turned on, close southern and western boundaries
#ifdef ALLOW_OBCS_SOUTH
         DO i=1-Olx,sNx+Olx
C Southern boundary
          J_obc = OB_Js(i,bi,bj)
          IF (J_obc.NE.0) THEN
#ifdef SEAICE_CGRID
           seaiceMaskU(i,J_obc,bi,bj)=   0.0 _d 0
           seaiceMaskV(i,J_obc,bi,bj)=   0.0 _d 0
#else
           UVM(i,J_obc,bi,bj)=0. _d 0
#endif /* SEAICE_CGRID */
          ENDIF
         ENDDO
#endif /* ALLOW_OBCS_SOUTH */
#ifdef ALLOW_OBCS_WEST
         DO j=1-Oly,sNy+Oly
C Western boundary
          I_obc=OB_Iw(j,bi,bj)
          IF (I_obc.NE.0) THEN
#ifdef SEAICE_CGRID
           seaiceMaskU(I_obc,j,bi,bj)=   0.0 _d 0
           seaiceMaskV(I_obc,j,bi,bj)=   0.0 _d 0
#else
           UVM(I_obc,j,bi,bj)=0. _d 0
#endif /* SEAICE_CGRID */
          ENDIF
         ENDDO
#endif /* ALLOW_OBCS_WEST */
        ENDIF
#endif /* ALLOW_OBCS */

#ifdef ALLOW_EXCH2
#ifndef SEAICE_CGRID
C--   Special stuff for cubed sphere: assume grid is rectangular and
C     set UV mask to zero except for Arctic and Antarctic cube faces.
        IF (useCubedSphereExchange) THEN
         myTile = W2_myTileList(bi,bj)
         IF ( exch2_myFace(myTile) .EQ. 1 .OR.
     &        exch2_myFace(myTile) .EQ. 2 .OR.
     &        exch2_myFace(myTile) .EQ. 4 .OR.
     &        exch2_myFace(myTile) .EQ. 5 ) THEN
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
            UVM(i,j,bi,bj)=0. _d 0
           ENDDO
          ENDDO
         ELSEIF ( exch2_isWedge(myTile) .EQ. 1 ) THEN
          i=1
          DO j=1-OLy,sNy+OLy
           UVM(i,j,bi,bj)=0. _d 0
          ENDDO
         ELSEIF ( exch2_isSedge(myTile) .EQ. 1 ) THEN
          j=1
          DO i=1-OLx,sNx+OLx
           UVM(i,j,bi,bj)=0. _d 0
          ENDDO
         ENDIF
        ENDIF
#endif /* SEAICE_CGRID */
#endif /* ALLOW_EXCH2 */

        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          TICE(i,j,bi,bj)=273.0 _d 0
#ifdef SEAICE_MULTICATEGORY
          DO k=1,MULTDIM
           TICES(i,j,k,bi,bj)=273.0 _d 0
          ENDDO
#endif /* SEAICE_MULTICATEGORY */
#ifndef SEAICE_CGRID
          AMASS      (i,j,bi,bj)=1000.0 _d 0
#else
          seaiceMassC(i,j,bi,bj)=1000.0 _d 0
          seaiceMassU(i,j,bi,bj)=1000.0 _d 0
          seaiceMassV(i,j,bi,bj)=1000.0 _d 0
#endif
         ENDDO
        ENDDO

       ENDDO
      ENDDO

C--   Update overlap regions
#ifdef SEAICE_CGRID
      CALL EXCH_UV_XY_RL(seaiceMaskU,seaiceMaskV,.FALSE.,myThid)
#else
      _EXCH_XY_RL(UVM, myThid)
#endif

C--   Now lets look at all these beasts
      IF ( debugLevel .GE. debLevB ) THEN
         myIter=0
         CALL PLOT_FIELD_XYRL( HEFFM   , 'Current HEFFM   ' ,
     &        myIter, myThid )
#ifdef SEAICE_CGRID
         CALL PLOT_FIELD_XYRL( seaiceMaskU, 'Current seaiceMaskU',
     &        myIter, myThid )
         CALL PLOT_FIELD_XYRL( seaiceMaskV, 'Current seaiceMaskV',
     &        myIter, myThid )
#else
         CALL PLOT_FIELD_XYRL( UVM     , 'Current UVM     ' ,
     &        myIter, myThid )
#endif
      ENDIF

C--   Set model variables to initial/restart conditions
      IF ( .NOT. ( startTime .EQ. baseTime .AND.  nIter0 .EQ. 0
     &     .AND. pickupSuff .EQ. ' ') ) THEN

         CALL SEAICE_READ_PICKUP ( myThid )

      ELSE

       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           TMIX(i,j,bi,bj)=TICE(i,j,bi,bj)
           HEFF(i,j,bi,bj)=SEAICE_initialHEFF*HEFFM(i,j,bi,bj)
           UICE(i,j,bi,bj)=ZERO
           VICE(i,j,bi,bj)=ZERO
          ENDDO
         ENDDO
        ENDDO
       ENDDO

C--   Read initial sea-ice thickness from file if available.
       IF ( HeffFile .NE. ' ' ) THEN
        CALL READ_FLD_XY_RL( HeffFile, ' ', HEFF, 0, myThid )
        _EXCH_XY_RL(HEFF,myThid)
        DO bj=myByLo(myThid),myByHi(myThid)
         DO bi=myBxLo(myThid),myBxHi(myThid)
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
            HEFF(i,j,bi,bj) = MAX(HEFF(i,j,bi,bj),ZERO)
           ENDDO
          ENDDO
         ENDDO
        ENDDO
       ENDIF

       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           IF(HEFF(i,j,bi,bj).GT.ZERO) AREA(i,j,bi,bj)=ONE
          ENDDO
         ENDDO
        ENDDO
       ENDDO

C--   Read initial sea-ice area from file if available.
       IF ( AreaFile .NE. ' ' ) THEN
        CALL READ_FLD_XY_RL( AreaFile, ' ', AREA, 0, myThid )
        _EXCH_XY_RL(AREA,myThid)
        DO bj=myByLo(myThid),myByHi(myThid)
         DO bi=myBxLo(myThid),myBxHi(myThid)
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
            AREA(i,j,bi,bj) = MAX(AREA(i,j,bi,bj),ZERO)
            AREA(i,j,bi,bj) = MIN(AREA(i,j,bi,bj),ONE)
            IF ( AREA(i,j,bi,bj) .LE. ZERO ) HEFF(i,j,bi,bj) = ZERO
            IF ( HEFF(i,j,bi,bj) .LE. ZERO ) AREA(i,j,bi,bj) = ZERO
           ENDDO
          ENDDO
         ENDDO
        ENDDO
       ENDIF

       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           HSNOW(i,j,bi,bj) = 0.2 _d 0 * AREA(i,j,bi,bj)
          ENDDO
         ENDDO
        ENDDO
       ENDDO

C--   Read initial snow thickness from file if available.
       IF ( HsnowFile .NE. ' ' ) THEN
        CALL READ_FLD_XY_RL( HsnowFile, ' ', HSNOW, 0, myThid )
        _EXCH_XY_RL(HSNOW,myThid)
        DO bj=myByLo(myThid),myByHi(myThid)
         DO bi=myBxLo(myThid),myBxHi(myThid)
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
            HSNOW(i,j,bi,bj) = MAX(HSNOW(i,j,bi,bj),ZERO)
           ENDDO
          ENDDO
         ENDDO
        ENDDO
       ENDIF

#ifdef SEAICE_VARIABLE_SALINITY
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           HSALT(i,j,bi,bj)=HEFF(i,j,bi,bj)*salt(i,j,kSurface,bi,bj)*
     &            SEAICE_rhoIce*SIsalFRAC
cif     &            ICE2WATR*rhoConstFresh*SIsalFRAC

          ENDDO
         ENDDO
        ENDDO
       ENDDO

C--   Read initial sea ice salinity from file if available.
       IF ( HsaltFile .NE. ' ' ) THEN
        CALL READ_FLD_XY_RL( HsaltFile, ' ', HSALT, 0, myThid )
        _EXCH_XY_RL(HSALT,myThid)
       ENDIF
#endif /* SEAICE_VARIABLE_SALINITY */

#ifdef SEAICE_AGE
C--   Read initial sea ice age from file if available.
       DO iTracer = 1, SEAICE_num
        IF ( IceAgeTrFile(iTracer) .NE. ' ' ) THEN
        CALL READ_FLD_XY_RL( IceAgeTrFile(iTracer), ' ',
     &   IceAgeTr(1-Olx,1-Oly,1,1,iTracer), 0, myThid )
        _EXCH_XY_RL(IceAgeTr(1-Olx,1-Oly,1,1,iTracer),myThid)
        ENDIF
       ENDDO
#endif /* SEAICE_AGE */

      ENDIF

C--   In case we use scheme with a large stencil that extends into overlap:
#ifdef ALLOW_OBCS
      IF ( useOBCS ) THEN
        DO bj=myByLo(myThid),myByHi(myThid)
         DO bi=myBxLo(myThid),myBxHi(myThid)
           CALL OBCS_COPY_TRACER( HEFF(1-Olx,1-Oly,bi,bj),
     I                            1, bi, bj, myThid )
           CALL OBCS_COPY_TRACER( AREA(1-Olx,1-Oly,bi,bj),
     I                            1, bi, bj, myThid )
           CALL OBCS_COPY_TRACER( HSNOW(1-Olx,1-Oly,bi,bj),
     I                            1, bi, bj, myThid )
#ifdef SEAICE_VARIABLE_SALINITY
           CALL OBCS_COPY_TRACER( HSALT(1-Olx,1-Oly,bi,bj),
     I                            1, bi, bj, myThid )
#endif
#ifdef SEAICE_AGE
           DO iTracer = 1, SEAICE_num
            CALL OBCS_COPY_TRACER( IceAgeTr(1-Olx,1-Oly,bi,bj,iTracer),
     I                             1, bi, bj, myThid )
           ENDDO
#endif
         ENDDO
        ENDDO
      ENDIF
#endif /* ALLOW_OBCS */

C---  Complete initialization
      PSTAR = SEAICE_strength
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          ZETA(i,j,bi,bj)   = HEFF(i,j,bi,bj)*(1.0 _d 11)
          ETA(i,j,bi,bj)    = ZETA(i,j,bi,bj)/SEAICE_eccen**2
          PRESS0(i,j,bi,bj) = PSTAR*HEFF(i,j,bi,bj)
     &         *EXP(-20.0 _d 0*(ONE-AREA(i,j,bi,bj)))
          ZMAX(I,J,bi,bj)   = SEAICE_zetaMaxFac*PRESS0(I,J,bi,bj)
          ZMIN(i,j,bi,bj)   = SEAICE_zetaMin
          PRESS0(i,j,bi,bj) = PRESS0(i,j,bi,bj)*HEFFM(i,j,bi,bj)
         ENDDO
        ENDDO
        IF ( useRealFreshWaterFlux .AND. .NOT.useThSIce ) THEN
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           sIceLoad(i,j,bi,bj) = HEFF(i,j,bi,bj)*SEAICE_rhoIce
     &                         + HSNOW(i,j,bi,bj)*SEAICE_rhoSnow

          ENDDO
         ENDDO
        ENDIF
       ENDDO
      ENDDO

      RETURN
      END
1	ifenty	1.56	C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_init_varia.F,v 1.55 2011/04/20 01:34:27 jmc Exp $
2	heimbach	1.1	C $Name: $
3
4			#include "SEAICE_OPTIONS.h"
5
6			CStartOfInterface
7			SUBROUTINE SEAICE_INIT_VARIA( myThid )
8	jmc	1.47	C ==========================================================
9	heimbach	1.1	C \| SUBROUTINE SEAICE_INIT_VARIA \|
10			C \| o Initialization of sea ice model. \|
11	jmc	1.47	C ==========================================================
12			C ==========================================================
13	heimbach	1.1	IMPLICIT NONE
14	jmc	1.23
15	heimbach	1.1	C === Global variables ===
16			#include "SIZE.h"
17			#include "EEPARAMS.h"
18			#include "PARAMS.h"
19			#include "GRID.h"
20	dimitri	1.21	#include "DYNVARS.h"
21	jmc	1.50	#include "FFIELDS.h"
22	heimbach	1.54	#include "SEAICE_SIZE.h"
23	jmc	1.50	#include "SEAICE_PARAMS.h"
24	heimbach	1.1	#include "SEAICE.h"
25	heimbach	1.54	#include "SEAICE_TRACER.h"
26	jmc	1.50	#include "SEAICE_TAVE.h"
27	heimbach	1.1	#ifdef ALLOW_EXCH2
28	jmc	1.38	# include "W2_EXCH2_SIZE.h"
29	dimitri	1.24	# include "W2_EXCH2_TOPOLOGY.h"
30			#endif
31	dimitri	1.11	#ifdef ALLOW_OBCS
32	dimitri	1.24	# include "OBCS_OPTIONS.h"
33			# include "OBCS.h"
34			#endif
35	heimbach	1.1
36			C === Routine arguments ===
37			C myThid - Thread no. that called this routine.
38			INTEGER myThid
39			CEndOfInterface
40	jmc	1.23
41	heimbach	1.1	C === Local variables ===
42			C i,j,k,bi,bj - Loop counters
43
44	jmc	1.55	INTEGER i, j, bi, bj
45	mlosch	1.8	_RL PSTAR
46	jmc	1.55	INTEGER myIter, kSurface
47			#ifdef SEAICE_CGRID
48	heimbach	1.1	_RS mask_uice
49	jmc	1.55	#endif
50	jmc	1.50	#ifdef SEAICE_MULTICATEGORY
51			INTEGER k
52			#endif
53	jmc	1.55	#ifdef SEAICE_AGE
54			INTEGER iTracer
55			#endif
56	dimitri	1.53	#ifdef ALLOW_OBCS
57			INTEGER I_obc, J_obc
58			#endif /* ALLOW_OBCS */
59	jmc	1.50	#ifdef ALLOW_EXCH2
60	dimitri	1.51	# ifndef SEAICE_CGRID
61	jmc	1.50	INTEGER myTile
62	dimitri	1.51	# endif
63	jmc	1.50	#endif
64	heimbach	1.1
65	jmc	1.23	IF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
66	mlosch	1.44	kSurface = Nr
67	jmc	1.23	ELSE
68	mlosch	1.44	kSurface = 1
69	jmc	1.23	ENDIF
70	dimitri	1.11
71	jmc	1.23	C-- Initialise all variables in common blocks:
72	heimbach	1.1	DO bj=myByLo(myThid),myByHi(myThid)
73			DO bi=myBxLo(myThid),myBxHi(myThid)
74	mlosch	1.42	DO j=1-OLy,sNy+OLy
75			DO i=1-OLx,sNx+OLx
76	mlosch	1.43	HEFF(i,j,bi,bj)=0. _d 0
77			AREA(i,j,bi,bj)=0. _d 0
78			UICE(i,j,bi,bj)=0. _d 0
79			VICE(i,j,bi,bj)=0. _d 0
80			C
81	mlosch	1.42	uIceNm1(i,j,bi,bj)=0. _d 0
82			vIceNm1(i,j,bi,bj)=0. _d 0
83			areaNm1(i,j,bi,bj)=0. _d 0
84			hEffNm1(i,j,bi,bj)=0. _d 0
85			ENDDO
86			ENDDO
87	jmc	1.23	#ifdef SEAICE_MULTICATEGORY
88			DO k=1,MULTDIM
89			DO j=1-OLy,sNy+OLy
90			DO i=1-OLx,sNx+OLx
91			TICES(i,j,k,bi,bj)=0. _d 0
92	heimbach	1.1	ENDDO
93			ENDDO
94			ENDDO
95	jmc	1.23	#endif
96			DO j=1-OLy,sNy+OLy
97			DO i=1-OLx,sNx+OLx
98			ETA (i,j,bi,bj) = 0. _d 0
99			ZETA(i,j,bi,bj) = 0. _d 0
100			DRAGS(i,j,bi,bj) = 0. _d 0
101			DRAGA(i,j,bi,bj) = 0. _d 0
102			FORCEX(i,j,bi,bj) = 0. _d 0
103			FORCEY(i,j,bi,bj) = 0. _d 0
104			UICEC(i,j,bi,bj) = 0. _d 0
105			VICEC(i,j,bi,bj) = 0. _d 0
106			#ifdef SEAICE_CGRID
107			seaiceMassC(i,j,bi,bj)=0. _d 0
108			seaiceMassU(i,j,bi,bj)=0. _d 0
109			seaiceMassV(i,j,bi,bj)=0. _d 0
110			stressDivergenceX(i,j,bi,bj) = 0. _d 0
111			stressDivergenceY(i,j,bi,bj) = 0. _d 0
112	heimbach	1.39	# ifdef SEAICE_ALLOW_EVP
113	jmc	1.23	seaice_sigma1 (i,j,bi,bj) = 0. _d 0
114			seaice_sigma2 (i,j,bi,bj) = 0. _d 0
115			seaice_sigma12(i,j,bi,bj) = 0. _d 0
116			# endif /* SEAICE_ALLOW_EVP */
117			#else /* SEAICE_CGRID */
118			AMASS(i,j,bi,bj) = 0. _d 0
119			DAIRN(i,j,bi,bj) = 0. _d 0
120			WINDX(i,j,bi,bj) = 0. _d 0
121			WINDY(i,j,bi,bj) = 0. _d 0
122	dimitri	1.25	GWATX(i,j,bi,bj) = 0. _d 0
123			GWATY(i,j,bi,bj) = 0. _d 0
124	jmc	1.23	#endif /* SEAICE_CGRID */
125			DWATN(i,j,bi,bj) = 0. _d 0
126			PRESS0(i,j,bi,bj) = 0. _d 0
127			FORCEX0(i,j,bi,bj)= 0. _d 0
128			FORCEY0(i,j,bi,bj)= 0. _d 0
129			ZMAX(i,j,bi,bj) = 0. _d 0
130			ZMIN(i,j,bi,bj) = 0. _d 0
131			HSNOW(i,j,bi,bj) = 0. _d 0
132	ifenty	1.56	#ifdef SEAICE_VARIABLE_SALINITY
133	jmc	1.23	HSALT(i,j,bi,bj) = 0. _d 0
134			#endif
135	dimitri	1.33	#ifdef SEAICE_AGE
136	heimbach	1.54	DO iTracer = 1, SEAICE_num
137			IceAgeTr(i,j,bi,bj,iTracer) = 0. _d 0
138			ENDDO
139	dimitri	1.33	#endif
140	jmc	1.23	YNEG (i,j,bi,bj) = 0. _d 0
141			RIVER(i,j,bi,bj) = 0. _d 0
142			TMIX(i,j,bi,bj) = 0. _d 0
143			TICE(i,j,bi,bj) = 0. _d 0
144			TAUX(i,j,bi,bj) = 0. _d 0
145			TAUY(i,j,bi,bj) = 0. _d 0
146			#ifdef ALLOW_SEAICE_COST_EXPORT
147			uHeffExportCell(i,j,bi,bj) = 0. _d 0
148			vHeffExportCell(i,j,bi,bj) = 0. _d 0
149	dimitri	1.18	#endif
150	heimbach	1.32	saltWtrIce(i,j,bi,bj) = 0. _d 0
151			frWtrIce(i,j,bi,bj) = 0. _d 0
152	heimbach	1.49	#ifdef ALLOW_MEAN_SFLUX_COST_CONTRIBUTION
153			frWtrAtm(i,j,bi,bj) = 0. _d 0
154			#endif
155	heimbach	1.1	ENDDO
156			ENDDO
157			ENDDO
158			ENDDO
159
160	jmc	1.50	#ifdef ALLOW_TIMEAVE
161			C Initialize averages to zero
162			DO bj = myByLo(myThid), myByHi(myThid)
163			DO bi = myBxLo(myThid), myBxHi(myThid)
164			CALL TIMEAVE_RESET( FUtave , 1, bi, bj, myThid )
165			CALL TIMEAVE_RESET( FVtave , 1, bi, bj, myThid )
166			CALL TIMEAVE_RESET( EmPmRtave, 1, bi, bj, myThid )
167			CALL TIMEAVE_RESET( QNETtave , 1, bi, bj, myThid )
168			CALL TIMEAVE_RESET( QSWtave , 1, bi, bj, myThid )
169			CALL TIMEAVE_RESET( UICEtave , 1, bi, bj, myThid )
170			CALL TIMEAVE_RESET( VICEtave , 1, bi, bj, myThid )
171			CALL TIMEAVE_RESET( HEFFtave , 1, bi, bj, myThid )
172			CALL TIMEAVE_RESET( AREAtave , 1, bi, bj, myThid )
173			SEAICE_timeAve(bi,bj) = ZERO
174			ENDDO
175			ENDDO
176			#endif /* ALLOW_TIMEAVE */
177
178	mlosch	1.40	C-- Initialize (variable) grid info. As long as we allow masking of
179			C-- velocities outside of ice covered areas (in seaice_dynsolver)
180			C-- we need to re-initialize seaiceMaskU/V here for TAF/TAMC
181	heimbach	1.1	DO bj=myByLo(myThid),myByHi(myThid)
182			DO bi=myBxLo(myThid),myBxHi(myThid)
183	jmc	1.23	DO j=1-OLy+1,sNy+OLy
184			DO i=1-OLx+1,sNx+OLx
185	heimbach	1.49	#ifdef SEAICE_CGRID
186	jmc	1.23	seaiceMaskU(i,j,bi,bj)= 0.0 _d 0
187			seaiceMaskV(i,j,bi,bj)= 0.0 _d 0
188			mask_uice=HEFFM(i,j,bi,bj)+HEFFM(i-1,j ,bi,bj)
189	heimbach	1.45	IF(mask_uice.GT.1.5 _d 0) seaiceMaskU(i,j,bi,bj)=1.0 _d 0
190	jmc	1.23	mask_uice=HEFFM(i,j,bi,bj)+HEFFM(i ,j-1,bi,bj)
191	heimbach	1.45	IF(mask_uice.GT.1.5 _d 0) seaiceMaskV(i,j,bi,bj)=1.0 _d 0
192	heimbach	1.49	#else
193			UVM(i,j,bi,bj)= 0.0 _d 0
194			#endif /* SEAICE_CGRID */
195	dimitri	1.22	ENDDO
196			ENDDO
197	mlosch	1.40	ENDDO
198			ENDDO
199	dimitri	1.11
200	mlosch	1.40	DO bj=myByLo(myThid),myByHi(myThid)
201			DO bi=myBxLo(myThid),myBxHi(myThid)
202	dimitri	1.11	#ifdef ALLOW_OBCS
203	dimitri	1.16	IF (useOBCS) THEN
204	dimitri	1.22	C-- If OBCS is turned on, close southern and western boundaries
205	dimitri	1.11	#ifdef ALLOW_OBCS_SOUTH
206	jmc	1.23	DO i=1-Olx,sNx+Olx
207	dimitri	1.11	C Southern boundary
208	dimitri	1.53	J_obc = OB_Js(i,bi,bj)
209			IF (J_obc.NE.0) THEN
210	dimitri	1.22	#ifdef SEAICE_CGRID
211	dimitri	1.53	seaiceMaskU(i,J_obc,bi,bj)= 0.0 _d 0
212			seaiceMaskV(i,J_obc,bi,bj)= 0.0 _d 0
213	dimitri	1.22	#else
214	dimitri	1.53	UVM(i,J_obc,bi,bj)=0. _d 0
215	dimitri	1.22	#endif /* SEAICE_CGRID */
216	dimitri	1.16	ENDIF
217			ENDDO
218	dimitri	1.22	#endif /* ALLOW_OBCS_SOUTH */
219	dimitri	1.11	#ifdef ALLOW_OBCS_WEST
220	jmc	1.23	DO j=1-Oly,sNy+Oly
221	dimitri	1.11	C Western boundary
222	dimitri	1.53	I_obc=OB_Iw(j,bi,bj)
223			IF (I_obc.NE.0) THEN
224	dimitri	1.22	#ifdef SEAICE_CGRID
225	dimitri	1.53	seaiceMaskU(I_obc,j,bi,bj)= 0.0 _d 0
226			seaiceMaskV(I_obc,j,bi,bj)= 0.0 _d 0
227	dimitri	1.22	#else
228	dimitri	1.53	UVM(I_obc,j,bi,bj)=0. _d 0
229	dimitri	1.22	#endif /* SEAICE_CGRID */
230	dimitri	1.16	ENDIF
231			ENDDO
232	dimitri	1.22	#endif /* ALLOW_OBCS_WEST */
233	dimitri	1.16	ENDIF
234	dimitri	1.11	#endif /* ALLOW_OBCS */
235
236	heimbach	1.1	#ifdef ALLOW_EXCH2
237	dimitri	1.11	#ifndef SEAICE_CGRID
238	heimbach	1.1	C-- Special stuff for cubed sphere: assume grid is rectangular and
239			C set UV mask to zero except for Arctic and Antarctic cube faces.
240			IF (useCubedSphereExchange) THEN
241	jmc	1.46	myTile = W2_myTileList(bi,bj)
242	heimbach	1.1	IF ( exch2_myFace(myTile) .EQ. 1 .OR.
243			& exch2_myFace(myTile) .EQ. 2 .OR.
244			& exch2_myFace(myTile) .EQ. 4 .OR.
245			& exch2_myFace(myTile) .EQ. 5 ) THEN
246	jmc	1.23	DO j=1-OLy,sNy+OLy
247			DO i=1-OLx,sNx+OLx
248	heimbach	1.10	UVM(i,j,bi,bj)=0. _d 0
249	heimbach	1.1	ENDDO
250			ENDDO
251			ELSEIF ( exch2_isWedge(myTile) .EQ. 1 ) THEN
252	jmc	1.23	i=1
253			DO j=1-OLy,sNy+OLy
254	heimbach	1.10	UVM(i,j,bi,bj)=0. _d 0
255	heimbach	1.1	ENDDO
256			ELSEIF ( exch2_isSedge(myTile) .EQ. 1 ) THEN
257	jmc	1.23	j=1
258			DO i=1-OLx,sNx+OLx
259	heimbach	1.10	UVM(i,j,bi,bj)=0. _d 0
260	heimbach	1.1	ENDDO
261			ENDIF
262			ENDIF
263	dimitri	1.11	#endif /* SEAICE_CGRID */
264	heimbach	1.1	#endif /* ALLOW_EXCH2 */
265
266			DO j=1-OLy,sNy+OLy
267			DO i=1-OLx,sNx+OLx
268	jmc	1.23	TICE(i,j,bi,bj)=273.0 _d 0
269	mlosch	1.2	#ifdef SEAICE_MULTICATEGORY
270	heimbach	1.1	DO k=1,MULTDIM
271	jmc	1.23	TICES(i,j,k,bi,bj)=273.0 _d 0
272	heimbach	1.1	ENDDO
273	mlosch	1.2	#endif /* SEAICE_MULTICATEGORY */
274	heimbach	1.1	#ifndef SEAICE_CGRID
275	jmc	1.23	AMASS (i,j,bi,bj)=1000.0 _d 0
276	heimbach	1.1	#else
277	jmc	1.23	seaiceMassC(i,j,bi,bj)=1000.0 _d 0
278			seaiceMassU(i,j,bi,bj)=1000.0 _d 0
279			seaiceMassV(i,j,bi,bj)=1000.0 _d 0
280	heimbach	1.1	#endif
281			ENDDO
282			ENDDO
283
284			ENDDO
285			ENDDO
286
287			C-- Update overlap regions
288			#ifdef SEAICE_CGRID
289			CALL EXCH_UV_XY_RL(seaiceMaskU,seaiceMaskV,.FALSE.,myThid)
290			#else
291	jmc	1.37	_EXCH_XY_RL(UVM, myThid)
292	heimbach	1.1	#endif
293
294			C-- Now lets look at all these beasts
295			IF ( debugLevel .GE. debLevB ) THEN
296			myIter=0
297			CALL PLOT_FIELD_XYRL( HEFFM , 'Current HEFFM ' ,
298			& myIter, myThid )
299			#ifdef SEAICE_CGRID
300			CALL PLOT_FIELD_XYRL( seaiceMaskU, 'Current seaiceMaskU',
301			& myIter, myThid )
302			CALL PLOT_FIELD_XYRL( seaiceMaskV, 'Current seaiceMaskV',
303			& myIter, myThid )
304			#else
305			CALL PLOT_FIELD_XYRL( UVM , 'Current UVM ' ,
306			& myIter, myThid )
307			#endif
308			ENDIF
309
310			C-- Set model variables to initial/restart conditions
311	mlosch	1.26	IF ( .NOT. ( startTime .EQ. baseTime .AND. nIter0 .EQ. 0
312			& .AND. pickupSuff .EQ. ' ') ) THEN
313	heimbach	1.1
314			CALL SEAICE_READ_PICKUP ( myThid )
315
316			ELSE
317
318			DO bj=myByLo(myThid),myByHi(myThid)
319			DO bi=myBxLo(myThid),myBxHi(myThid)
320			DO j=1-OLy,sNy+OLy
321			DO i=1-OLx,sNx+OLx
322	jmc	1.23	TMIX(i,j,bi,bj)=TICE(i,j,bi,bj)
323	mlosch	1.43	HEFF(i,j,bi,bj)=SEAICE_initialHEFF*HEFFM(i,j,bi,bj)
324			UICE(i,j,bi,bj)=ZERO
325			VICE(i,j,bi,bj)=ZERO
326	heimbach	1.1	ENDDO
327			ENDDO
328			ENDDO
329			ENDDO
330
331			C-- Read initial sea-ice thickness from file if available.
332			IF ( HeffFile .NE. ' ' ) THEN
333	mlosch	1.43	CALL READ_FLD_XY_RL( HeffFile, ' ', HEFF, 0, myThid )
334			_EXCH_XY_RL(HEFF,myThid)
335	heimbach	1.1	DO bj=myByLo(myThid),myByHi(myThid)
336			DO bi=myBxLo(myThid),myBxHi(myThid)
337			DO j=1-OLy,sNy+OLy
338			DO i=1-OLx,sNx+OLx
339	mlosch	1.43	HEFF(i,j,bi,bj) = MAX(HEFF(i,j,bi,bj),ZERO)
340	heimbach	1.1	ENDDO
341			ENDDO
342			ENDDO
343			ENDDO
344			ENDIF
345
346			DO bj=myByLo(myThid),myByHi(myThid)
347			DO bi=myBxLo(myThid),myBxHi(myThid)
348			DO j=1-OLy,sNy+OLy
349			DO i=1-OLx,sNx+OLx
350	mlosch	1.43	IF(HEFF(i,j,bi,bj).GT.ZERO) AREA(i,j,bi,bj)=ONE
351	heimbach	1.1	ENDDO
352			ENDDO
353			ENDDO
354			ENDDO
355	jmc	1.23
356	dimitri	1.27	C-- Read initial sea-ice area from file if available.
357	mlosch	1.7	IF ( AreaFile .NE. ' ' ) THEN
358	mlosch	1.43	CALL READ_FLD_XY_RL( AreaFile, ' ', AREA, 0, myThid )
359			_EXCH_XY_RL(AREA,myThid)
360	mlosch	1.7	DO bj=myByLo(myThid),myByHi(myThid)
361			DO bi=myBxLo(myThid),myBxHi(myThid)
362			DO j=1-OLy,sNy+OLy
363			DO i=1-OLx,sNx+OLx
364	mlosch	1.43	AREA(i,j,bi,bj) = MAX(AREA(i,j,bi,bj),ZERO)
365			AREA(i,j,bi,bj) = MIN(AREA(i,j,bi,bj),ONE)
366			IF ( AREA(i,j,bi,bj) .LE. ZERO ) HEFF(i,j,bi,bj) = ZERO
367			IF ( HEFF(i,j,bi,bj) .LE. ZERO ) AREA(i,j,bi,bj) = ZERO
368	mlosch	1.7	ENDDO
369			ENDDO
370			ENDDO
371			ENDDO
372			ENDIF
373
374			DO bj=myByLo(myThid),myByHi(myThid)
375			DO bi=myBxLo(myThid),myBxHi(myThid)
376			DO j=1-OLy,sNy+OLy
377			DO i=1-OLx,sNx+OLx
378	jmc	1.47	HSNOW(i,j,bi,bj) = 0.2 _d 0 * AREA(i,j,bi,bj)
379	mlosch	1.7	ENDDO
380			ENDDO
381			ENDDO
382			ENDDO
383	dimitri	1.9
384			C-- Read initial snow thickness from file if available.
385			IF ( HsnowFile .NE. ' ' ) THEN
386	mlosch	1.43	CALL READ_FLD_XY_RL( HsnowFile, ' ', HSNOW, 0, myThid )
387			_EXCH_XY_RL(HSNOW,myThid)
388	dimitri	1.9	DO bj=myByLo(myThid),myByHi(myThid)
389			DO bi=myBxLo(myThid),myBxHi(myThid)
390			DO j=1-OLy,sNy+OLy
391			DO i=1-OLx,sNx+OLx
392	mlosch	1.43	HSNOW(i,j,bi,bj) = MAX(HSNOW(i,j,bi,bj),ZERO)
393	dimitri	1.9	ENDDO
394			ENDDO
395			ENDDO
396			ENDDO
397			ENDIF
398	dimitri	1.18
399	ifenty	1.56	#ifdef SEAICE_VARIABLE_SALINITY
400	dimitri	1.20	DO bj=myByLo(myThid),myByHi(myThid)
401			DO bi=myBxLo(myThid),myBxHi(myThid)
402	jmc	1.23	DO j=1-OLy,sNy+OLy
403			DO i=1-OLx,sNx+OLx
404	mlosch	1.44	HSALT(i,j,bi,bj)=HEFF(i,j,bi,bj)salt(i,j,kSurface,bi,bj)
405	ifenty	1.56	& SEAICE_rhoIce*SIsalFRAC
406			cif & ICE2WATRrhoConstFreshSIsalFRAC
407
408	dimitri	1.20	ENDDO
409			ENDDO
410			ENDDO
411			ENDDO
412
413	dimitri	1.18	C-- Read initial sea ice salinity from file if available.
414			IF ( HsaltFile .NE. ' ' ) THEN
415	mlosch	1.43	CALL READ_FLD_XY_RL( HsaltFile, ' ', HSALT, 0, myThid )
416			_EXCH_XY_RL(HSALT,myThid)
417	dimitri	1.18	ENDIF
418	ifenty	1.56	#endif /* SEAICE_VARIABLE_SALINITY */
419	jmc	1.23
420	dimitri	1.33	#ifdef SEAICE_AGE
421			C-- Read initial sea ice age from file if available.
422	heimbach	1.54	DO iTracer = 1, SEAICE_num
423			IF ( IceAgeTrFile(iTracer) .NE. ' ' ) THEN
424	jmc	1.55	CALL READ_FLD_XY_RL( IceAgeTrFile(iTracer), ' ',
425	heimbach	1.54	& IceAgeTr(1-Olx,1-Oly,1,1,iTracer), 0, myThid )
426			_EXCH_XY_RL(IceAgeTr(1-Olx,1-Oly,1,1,iTracer),myThid)
427			ENDIF
428			ENDDO
429	dimitri	1.33	#endif /* SEAICE_AGE */
430
431	heimbach	1.1	ENDIF
432
433	jmc	1.52	C-- In case we use scheme with a large stencil that extends into overlap:
434			#ifdef ALLOW_OBCS
435			IF ( useOBCS ) THEN
436			DO bj=myByLo(myThid),myByHi(myThid)
437			DO bi=myBxLo(myThid),myBxHi(myThid)
438			CALL OBCS_COPY_TRACER( HEFF(1-Olx,1-Oly,bi,bj),
439			I 1, bi, bj, myThid )
440			CALL OBCS_COPY_TRACER( AREA(1-Olx,1-Oly,bi,bj),
441			I 1, bi, bj, myThid )
442			CALL OBCS_COPY_TRACER( HSNOW(1-Olx,1-Oly,bi,bj),
443			I 1, bi, bj, myThid )
444	ifenty	1.56	#ifdef SEAICE_VARIABLE_SALINITY
445	jmc	1.52	CALL OBCS_COPY_TRACER( HSALT(1-Olx,1-Oly,bi,bj),
446			I 1, bi, bj, myThid )
447			#endif
448			#ifdef SEAICE_AGE
449	heimbach	1.54	DO iTracer = 1, SEAICE_num
450			CALL OBCS_COPY_TRACER( IceAgeTr(1-Olx,1-Oly,bi,bj,iTracer),
451			I 1, bi, bj, myThid )
452			ENDDO
453	jmc	1.52	#endif
454			ENDDO
455			ENDDO
456			ENDIF
457			#endif /* ALLOW_OBCS */
458
459	heimbach	1.1	C--- Complete initialization
460	mlosch	1.8	PSTAR = SEAICE_strength
461	heimbach	1.1	DO bj=myByLo(myThid),myByHi(myThid)
462			DO bi=myBxLo(myThid),myBxHi(myThid)
463			DO j=1-OLy,sNy+OLy
464			DO i=1-OLx,sNx+OLx
465	mlosch	1.48	ZETA(i,j,bi,bj) = HEFF(i,j,bi,bj)*(1.0 _d 11)
466			ETA(i,j,bi,bj) = ZETA(i,j,bi,bj)/SEAICE_eccen**2
467			PRESS0(i,j,bi,bj) = PSTAR*HEFF(i,j,bi,bj)
468	mlosch	1.43	& EXP(-20.0 _d 0(ONE-AREA(i,j,bi,bj)))
469	mlosch	1.48	ZMAX(I,J,bi,bj) = SEAICE_zetaMaxFac*PRESS0(I,J,bi,bj)
470			ZMIN(i,j,bi,bj) = SEAICE_zetaMin
471			PRESS0(i,j,bi,bj) = PRESS0(i,j,bi,bj)*HEFFM(i,j,bi,bj)
472	heimbach	1.1	ENDDO
473			ENDDO
474			IF ( useRealFreshWaterFlux .AND. .NOT.useThSIce ) THEN
475			DO j=1-OLy,sNy+OLy
476			DO i=1-OLx,sNx+OLx
477	mlosch	1.43	sIceLoad(i,j,bi,bj) = HEFF(i,j,bi,bj)*SEAICE_rhoIce
478	mlosch	1.36	& + HSNOW(i,j,bi,bj)*SEAICE_rhoSnow
479	jmc	1.23
480	heimbach	1.1	ENDDO
481			ENDDO
482			ENDIF
483			ENDDO
484			ENDDO
485
486			RETURN
487			END