C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/arctic40km/code/Attic/seaice_init.F,v 1.2 2006/05/06 13:46:39 dimitri dead $
C $Name:  $

#include "SEAICE_OPTIONS.h"

CStartOfInterface
      SUBROUTINE SEAICE_INIT( myThid )
C     /==========================================================\
C     | SUBROUTINE SEAICE_INIT                                   |
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 "SEAICE.h"
#include "SEAICE_GRID.h"
#include "SEAICE_DIAGS.h"
#include "SEAICE_PARAMS.h"
#ifdef ALLOW_EXCH2
#include "W2_EXCH2_TOPOLOGY.h"
#include "W2_EXCH2_PARAMS.h"
#endif /* ALLOW_EXCH2 */

C     === Routine arguments ===
C     myThid - Thread no. that called this routine.
      INTEGER myThid
CEndOfInterface
 
#ifdef ALLOW_SEAICE
C     === Local variables ===
C     i,j,k,bi,bj - Loop counters

      INTEGER i, j, k, bi, bj
      _RS  mask_uice
      INTEGER myIter, myTile
      INTEGER iG, jG
      _RL lat, dlat, dlon, xG0, yG0
      _RL xGloc(1-Olx:sNx+Olx+1,1-Oly:sNy+Oly+1)
      _RL yGloc(1-Olx:sNx+Olx+1,1-Oly:sNy+Oly+1)
      INTEGER iGl,jGl
      iGl(I,bi) = 1+mod(myXGlobalLo-1+(bi-1)*sNx+I+Olx*Nx-1,Nx)
      jGl(J,bj) = 1+mod(myYGlobalLo-1+(bj-1)*sNy+J+Oly*Ny-1,Ny)

#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)
        DO k=1,Nr
         SEAICE_TimeAve(k,bi,bj)=ZERO
        ENDDO
       ENDDO
      ENDDO
#endif /* ALLOW_TIMEAVE */

cph(
cph   make sure TAF sees proper initialisation
cph   to avoid partial recomputation issues
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
c
        DO K=1,3
         DO J=1-OLy,sNy+OLy
          DO I=1-OLx,sNx+OLx
           HEFF(I,J,k,bi,bj)=ZERO
           AREA(I,J,k,bi,bj)=ZERO
           UICE(I,J,k,bi,bj)=ZERO
           VICE(I,J,k,bi,bj)=ZERO
          ENDDO
         ENDDO
        ENDDO
c
        DO J=1-OLy,sNy+OLy
         DO I=1-OLx,sNx+OLx
          HSNOW(I,J,bi,bj)=ZERO
          ZETA(I,J,bi,bj)=ZERO
         ENDDO
        ENDDO
c
       ENDDO
      ENDDO
cph)

C_JZ- The following is for a rotated model grid (for the UW regional arctic model)
C     For each tile ...
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)

C--     "Global" index (place holder)
        jG = myYGlobalLo + (bj-1)*sNy
        iG = myXGlobalLo + (bi-1)*sNx

C--   First find coordinate of tile corner (meaning outer corner of halo)
        xG0 = thetaMin
C       Find the X-coordinate of the outer grid-line of the "real" tile
        DO i=1, iG-1
         xG0 = xG0 + delX(i)
        ENDDO
C       Back-step to the outer grid-line of the "halo" region
        DO i=1, Olx
         xG0 = xG0 - delX( 1+mod(Olx*Nx-1+iG-i,Nx) )
        ENDDO
C       Find the Y-coordinate of the outer grid-line of the "real" tile
        yG0 = phiMin
        DO j=1, jG-1
         yG0 = yG0 + delY(j)
        ENDDO
C       Back-step to the outer grid-line of the "halo" region
        DO j=1, Oly
         yG0 = yG0 - delY( 1+mod(Oly*Ny-1+jG-j,Ny) )
        ENDDO

C--     Calculate coordinates of cell corners for N+1 grid-lines
        DO J=1-Oly,sNy+Oly +1
         xGloc(1-Olx,J) = xG0
         DO I=1-Olx,sNx+Olx
c         xGloc(I+1,J) = xGloc(I,J) + delX(1+mod(Nx-1+iG-1+i,Nx))
          xGloc(I+1,J) = xGloc(I,J) + delX( iGl(I,bi) )
         ENDDO
        ENDDO
        DO I=1-Olx,sNx+Olx +1
         yGloc(I,1-Oly) = yG0
         DO J=1-Oly,sNy+Oly
c         yGloc(I,J+1) = yGloc(I,J) + delY(1+mod(Ny-1+jG-1+j,Ny))
          yGloc(I,J+1) = yGloc(I,J) + delY( jGl(J,bj) )
         ENDDO
        ENDDO
        DO J=1-Oly,sNy+Oly
         DO I=1-Olx,sNx+Olx
          CSUICE(i,j,bi,bj) =cos(yGloc(I,J)*deg2rad)
          SINEICE(i,j,bi,bj) =sin(yGloc(I,J)*deg2rad)
          TNGICE(i,j,bi,bj)=SINEICE(i,j,bi,bj)/CSUICE(i,j,bi,bj)
         ENDDO
        ENDDO

C--     Calculate [xC,yC], coordinates of cell centers by averaging
        DO J=1-Oly,sNy+Oly
         DO I=1-Olx,sNx+Olx
          CSTICE(i,j,bi,bj) =cos(0.25*(
     &     yGloc(I,J)+yGloc(I+1,J)+yGloc(I,J+1)+yGloc(I+1,J+1))*deg2rad)
          SINEICE(i,j,bi,bj) =sin(0.25*(
     &     yGloc(I,J)+yGloc(I+1,J)+yGloc(I,J+1)+yGloc(I+1,J+1))*deg2rad)
          TNGTICE(i,j,bi,bj)=SINEICE(i,j,bi,bj)/CSTICE(i,j,bi,bj)
         ENDDO
        ENDDO

       ENDDO ! bi
      ENDDO ! bj

C--   Initialize grid info
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO J=1-OLy,sNy+OLy
         DO I=1-OLx,sNx+OLx
          SINEICE(i,j,bi,bj)=sin(yG(I,J,bi,bj)*deg2rad)
          DXTICE(i,j,bi,bj)=dxF(i,j,bi,bj)/CSTICE(i,j,bi,bj)
          DXUICE(i,j,bi,bj)=dxV(i,j,bi,bj)/CSUICE(i,j,bi,bj)
          DYTICE(i,j,bi,bj)=dyF(i,j,bi,bj)
          DYUICE(i,j,bi,bj)=dyU(i,j,bi,bj)
         ENDDO
        ENDDO
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          HEFFM(i,j,bi,bj)=ONE
          IF (_hFacC(i,j,1,bi,bj).eq.0.) HEFFM(i,j,bi,bj)=ZERO
         ENDDO
        ENDDO
        DO J=2-OLy,sNy+OLy
         DO I=2-OLx,sNx+OLx
          UVM(i,j,bi,bj)=ZERO
          mask_uice=HEFFM(I,J,  bi,bj)+HEFFM(I-1,J-1,bi,bj)
     &             +HEFFM(I,J-1,bi,bj)+HEFFM(I-1,J,  bi,bj)
          IF(mask_uice.GT.3.5) UVM(I,J,bi,bj)=ONE
         ENDDO
        ENDDO

#ifdef ALLOW_EXCH2
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
         DO J=1-OLy,sNy+OLy
          DO I=1-OLx,sNx+OLx
           CSTICE(i,j,bi,bj) = ONE
           CSUICE(i,j,bi,bj) = ONE
           TNGTICE(i,j,bi,bj)= ZERO
           TNGICE(i,j,bi,bj) = ZERO
           DXTICE(i,j,bi,bj) = dxF(i,j,bi,bj)
           DXUICE(i,j,bi,bj) = dxV(i,j,bi,bj)
          ENDDO
         ENDDO
         myTile = W2_myTileList(bi)
         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)=ZERO
           ENDDO
          ENDDO
         ELSEIF ( exch2_isWedge(myTile) .EQ. 1 ) THEN
          I=1
          DO J=1-OLy,sNy+OLy
           UVM(i,j,bi,bj)=ZERO
          ENDDO
         ELSEIF ( exch2_isSedge(myTile) .EQ. 1 ) THEN
          J=1
          DO I=1-OLx,sNx+OLx
           UVM(i,j,bi,bj)=ZERO
          ENDDO
         ENDIF
        ENDIF
#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_MULTILEVEL
          DO k=1,MULTDIM
           TICES(I,J,k,bi,bj)=273.0 _d 0
          ENDDO
#endif
          UICEC(I,J,bi,bj)=ZERO
          VICEC(I,J,bi,bj)=ZERO
          AMASS(I,J,bi,bj)=1000.0 _d 0
         ENDDO
        ENDDO
       ENDDO
      ENDDO

C--   Update overlap regions
      _EXCH_XY_R8(UVM, myThid)

      myIter=0
      CALL PLOT_FIELD_XYRL( CSTICE  , 'Current CSTICE  ' ,
     &                        myIter, myThid )
      CALL PLOT_FIELD_XYRL( CSUICE  , 'Current CSUICE  ' ,
     &                        myIter, myThid )
      CALL PLOT_FIELD_XYRL( TNGTICE  , 'Current TNGTICE  ' ,
     &                        myIter, myThid )
      CALL PLOT_FIELD_XYRL( TNGICE  , 'Current TNGICE  ' ,
     &                        myIter, myThid )
      CALL PLOT_FIELD_XYRL( SINEICE  , 'Current SINEICE  ' ,
     &                        myIter, myThid )
      CALL PLOT_FIELD_XYRL( DXTICE  , 'Current DXTICE  ' ,
     &                        myIter, myThid )
      CALL PLOT_FIELD_XYRL( DXUICE  , 'Current DXUICE  ' ,
     &                        myIter, myThid )
      CALL PLOT_FIELD_XYRL( DYTICE  , 'Current DYTICE  ' ,
     &                        myIter, myThid )
      CALL PLOT_FIELD_XYRL( DYUICE  , 'Current DYUICE  ' ,
     &                        myIter, myThid )
      CALL PLOT_FIELD_XYRL( HEFFM  , 'Current HEFFM  ' ,
     &                        myIter, myThid )
      CALL PLOT_FIELD_XYRL( UVM  , 'Current UVM  ' ,
     &                        myIter, myThid )

C--   Set model variables to initial/restart conditions
      IF ( nIter0 .NE. 0 ) 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
           HSNOW(I,J,bi,bj)=0.2*HEFFM(i,j,bi,bj)
           YNEG(I,J,bi,bj)=ZERO
           TMIX(I,J,bi,bj)=TICE(I,J,bi,bj)
           DO k=1,3
            HEFF(I,J,k,bi,bj)=SEAICE_initialHEFF*HEFFM(i,j,bi,bj)
            UICE(I,J,k,bi,bj)=ZERO
            VICE(I,J,k,bi,bj)=ZERO
           ENDDO
          ENDDO
         ENDDO
        ENDDO
       ENDDO

C--   Read initial sea-ice thickness from file if available.
       IF ( HeffFile .NE. ' ' ) THEN
        _BEGIN_MASTER( myThid )
        CALL READ_FLD_XY_RL( HeffFile, ' ', ZETA, 0, myThid )
        _END_MASTER(myThid)
        _EXCH_XY_R8(ZETA,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
            DO k=1,3
             HEFF(I,J,k,bi,bj) = MAX(ZETA(i,j,bi,bj),ZERO)
            ENDDO
           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
           DO k=1,3
            IF(HEFF(I,J,k,bi,bj).GT.ZERO) AREA(I,J,k,bi,bj)=ONE
           ENDDO
          ENDDO
         ENDDO
        ENDDO
       ENDDO
         
      ENDIF

C---  Complete initialization
      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,1,bi,bj)*(1.0 _d 11)
          ETA(I,J,bi,bj)=ZETA(I,J,bi,bj)/4.0 _d 0
         ENDDO
        ENDDO
       ENDDO
      ENDDO

#endif /* ALLOW_SEAICE */

      RETURN
      END