#include "SEAICE_OPTIONS.h"

CBOP
C     !ROUTINE: SEAICE_DIAGNOSTICS_STATE
C     !INTERFACE:
      SUBROUTINE SEAICE_DIAGNOSTICS_STATE(
     I                      myTime, myIter, myThid )
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | S/R  SEAICE_DIAGNOSTICS_STATE
C     | o fill-in diagnostics array for SEAICE state variables
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE

C     == Global variables ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "SEAICE_SIZE.h"
#include "SEAICE_PARAMS.h"
#include "SEAICE.h"
#include "SEAICE_TRACER.h"

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine Arguments ==
C     myTime  :: time counter for this thread
C     myIter  :: iteration counter for this thread
C     bi,bj   :: tile indices
C     myThid  :: thread number for this instance of the routine.
      _RL  myTime
      INTEGER myIter
      INTEGER myThid
CEOP

#ifdef ALLOW_DIAGNOSTICS
C     == Local variables ==
      INTEGER i,j,bi,bj
      _RL tIce (1-oLx:sNx+oLx,1-oLy:sNy+oLy,nSx,nSy)
#ifdef SEAICE_CGRID
      _RL sig1(1-oLx:sNx+oLx,1-oLy:sNy+oLy)
      _RL sig2(1-oLx:sNx+oLx,1-oLy:sNy+oLy)
      _RL sig11(1-oLx:sNx+oLx,1-oLy:sNy+oLy)
      _RL sig22(1-oLx:sNx+oLx,1-oLy:sNy+oLy)
      _RL sig12(1-oLx:sNx+oLx,1-oLy:sNy+oLy)
      _RL sigp, sigm, sigTmp, recip_prs
#endif
      INTEGER k
      _RL recip_multDim
      _RL tmp
#ifdef ALLOW_SITRACER
      INTEGER iTracer
      CHARACTER*8   diagName
#endif

      LOGICAL  DIAGNOSTICS_IS_ON
      EXTERNAL DIAGNOSTICS_IS_ON

      IF ( useDiagnostics ) THEN

       CALL DIAGNOSTICS_FILL( AREA, 'SIarea  ', 0, 1, 0, 1,1, myThid )
       CALL DIAGNOSTICS_FILL( HEFF, 'SIheff  ', 0, 1, 0, 1,1, myThid )
       CALL DIAGNOSTICS_FILL( UICE, 'SIuice  ', 0, 1, 0, 1,1, myThid )
       CALL DIAGNOSTICS_FILL( VICE, 'SIvice  ', 0, 1, 0, 1,1, myThid )

       IF ( DIAGNOSTICS_IS_ON('SItices ',myThid) ) THEN
        DO bj = myByLo(myThid), myByHi(myThid)
         DO bi = myBxLo(myThid), myBxHi(myThid)
C     use tIce as a temporary field, as it is done in seaice_growth
          DO j=1,sNy
           DO i=1,sNx
            tIce(I,J,bi,bj) = 0. _d 0
           ENDDO
          ENDDO
C     division by zero is not possible
          recip_multDim = 1. _d 0/SEAICE_multDim
          DO k=1,SEAICE_multDim
           DO j=1,sNy
            DO i=1,sNx
             tmp = 1. _d 0
#ifdef SEAICE_ITD
             IF (AREA(I,J,bi,bj) .GT. ZERO) 
     &            tmp=AREAITD(I,J,K,bi,bj)/AREA(I,J,bi,bj)
#endif /* SEAICE_ITD */
             tIce(I,J,bi,bj) = tIce(I,J,bi,bj)
     &            + TICES(I,J,K,bi,bj)*tmp*recip_multDim
            ENDDO
           ENDDO
          ENDDO
         ENDDO
        ENDDO
        CALL DIAGNOSTICS_FRACT_FILL(
     I                   tIce, AREA, 1. _d 0, 1, 'SItices ',
     I                   0, 1, 0, 1, 1, myThid )
       ENDIF

       CALL DIAGNOSTICS_FILL(HSNOW  ,'SIhsnow ',0,1 ,0,1,1,myThid)

#ifdef SEAICE_ITD
       CALL DIAGNOSTICS_FILL(HEFFITD ,'SIheffN ',0,nITD ,0,1,1,myThid)
       CALL DIAGNOSTICS_FILL(AREAITD ,'SIareaN ',0,nITD ,0,1,1,myThid)
       CALL DIAGNOSTICS_FILL(HSNOWITD,'SIhsnowN',0,nITD ,0,1,1,myThid)
#endif

#ifdef ALLOW_SITRACER
      DO iTracer = 1, SItrNumInUse
       WRITE(diagName,'(A4,I2.2,A2)') 'SItr',iTracer,'  '
       if (SItrMate(iTracer).EQ.'HEFF') then
        CALL DIAGNOSTICS_FRACT_FILL(
     I       SItracer(1-OLx,1-OLy,1,1,iTracer), HEFF,
     I       ONE, 1, diagName, 0, 1, 0, 1, 1, myThid )
       else
        CALL DIAGNOSTICS_FRACT_FILL(
     I       SItracer(1-OLx,1-OLy,1,1,iTracer), AREA,
     I       ONE, 1, diagName, 0, 1, 0, 1, 1, myThid )
       endif
      ENDDO
#endif
#ifdef SEAICE_VARIABLE_SALINITY
       CALL DIAGNOSTICS_FILL(HSALT  ,'SIhsalt ',0,1 ,0,1,1,myThid)
#endif

       CALL DIAGNOSTICS_FILL(zeta   ,'SIzeta  ',0,1 ,0,1,1,myThid)
       CALL DIAGNOSTICS_FILL(eta    ,'SIeta   ',0,1 ,0,1,1,myThid)
       CALL DIAGNOSTICS_FILL(press  ,'SIpress ',0,1 ,0,1,1,myThid)
       CALL DIAGNOSTICS_FILL(deltaC ,'SIdelta ',0,1 ,0,1,1,myThid)
#ifdef SEAICE_CGRID
       IF ( DIAGNOSTICS_IS_ON('SItensil',myThid) ) THEN
        DO bj = myByLo(myThid), myByHi(myThid)
         DO bi = myBxLo(myThid), myBxHi(myThid)
C     use tIce as a temporary field
          DO j=1,sNy
           DO i=1,sNx
            IF ( tensileStrFac(I,J,bi,bj) .EQ. 1. _d 0) THEN
C     This special case of tensile strength equal to compressive strength
C     is not very physical and should actually not happen but you never know;
C     in this case, press = P-T = P*(1-k) = 0. and we have to use press0 to
C     get something
             tIce(I,J,bi,bj) = press0(I,J,bi,bj)
            ELSE
C     This is more complicated than you think because press = P-T = P*(1-k),
C     but we are looking for T = k*P = k*press/(1-k)
             tIce(I,J,bi,bj) = tensileStrFac(I,J,bi,bj)
     &            *press(I,J,bi,bj)/(1. _d 0 - tensileStrFac(I,J,bi,bj))
            ENDIF
           ENDDO
          ENDDO
         ENDDO
        ENDDO
        CALL DIAGNOSTICS_FILL(tIce,'SItensil',0,1 ,0,1,1,myThid)
       ENDIF
       IF ( DIAGNOSTICS_IS_ON('SIsig1  ',myThid) .OR.
     &      DIAGNOSTICS_IS_ON('SIsig2  ',myThid) .OR.
     &      DIAGNOSTICS_IS_ON('SIshear ',myThid) ) THEN
#ifdef SEAICE_ALLOW_EVP
        IF ( SEAICEuseEVP ) THEN
C     for EVP compute principle stress components from recent
C     stress state and normalize with latest
C     PRESS = PRESS(n-1), n = number of sub-cycling steps
         DO bj = myByLo(myThid), myByHi(myThid)
          DO bi = myBxLo(myThid), myBxHi(myThid)
           DO j=1,sNy
            DO i=1,sNx
             sigp = seaice_sigma1(I,J,bi,bj)
             sigm = seaice_sigma2(I,J,bi,bj)
             sig12(I,J) = 0.25 _d 0 *
     &            ( seaice_sigma12(I,  J,  bi,bj)
     &            + seaice_sigma12(I+1,J,  bi,bj)
     &            + seaice_sigma12(I+1,J+1,bi,bj)
     &            + seaice_sigma12(I  ,J+1,bi,bj) )
             sigTmp = SQRT( sigm*sigm + 4. _d 0*sig12(I,J)*sig12(I,J) )
             recip_prs = 0. _d 0
             IF ( press0(I,J,bi,bj) .GT. 1. _d -13 )
     &            recip_prs = 1./press0(I,J,bi,bj)
             sig1(I,J) = 0.5*(sigp + sigTmp)*recip_prs
             sig2(I,J) = 0.5*(sigp - sigTmp)*recip_prs
            ENDDO
           ENDDO
           CALL DIAGNOSTICS_FILL(sig1,'SIsig1  ',0,1,2,bi,bj,myThid)
           CALL DIAGNOSTICS_FILL(sig2,'SIsig2  ',0,1,2,bi,bj,myThid)
          ENDDO
         ENDDO
        ELSE
#else
        IF ( .TRUE. ) THEN
#endif /* SEAICE_ALLOW_EVP */
C     recompute strainrates from up-to-date velocities
         CALL SEAICE_CALC_STRAINRATES(
     I        uIce, vIce,
     O        e11, e22, e12,
     I        0, myTime, myIter, myThid )
CML      CALL SEAICE_CALC_VISCOSITIES(
CML     I     e11, e22, e12, zMin, zMax, HEFFM, press0, tensileStrFac,
CML     O     eta, etaZ, zeta, zetaZ, press, deltaC,
CML     I     0, myTime, myIter, myThid )
C     but use old viscosities and pressure for the
C     principle stress components
         DO bj = myByLo(myThid), myByHi(myThid)
          DO bi = myBxLo(myThid), myBxHi(myThid)
           CALL SEAICE_CALC_STRESS(
     I          e11, e22, e12, press, zeta, eta, etaZ,
     O          sig11, sig22, sig12,
     I          bi, bj, myTime, myIter, myThid )
           DO j=1,sNy
            DO i=1,sNx
             sigp   = sig11(I,J) + sig22(I,J)
             sigm   = sig11(I,J) - sig22(I,J)
C     This should be the way of computing sig12 at C-points, 
C            sigTmp = 0.25 _d 0 *
C    &            ( sig12(I  ,J  ) + sig12(I+1,J  )
C    &            + sig12(I  ,J+1) + sig12(I+1,J+1) )
C     but sig12 = 2*etaZ*e12, and because of strong gradients in eta, 
C     etaZ can be very large for a cell with small eta and the straightforward
C     way of averaging mixes large etaZ with small press0, so we have to do it
C     in different way to get meaningfull sig12C (=sigTmp):
             sigTmp = 2.*eta(I,J,bi,bj) * 0.25 _d 0 *
     &            (e12(I,J,bi,bj) + e12(I+1,J,bi,bj)
     &            +e12(I,J+1,bi,bj)+e12(I+1,J+1,bi,bj))
             sigTmp = SQRT( sigm*sigm + 4.*sigTmp*sigTmp )
             recip_prs = 0. _d 0
             IF ( press0(I,J,bi,bj) .GT. 1. _d -13 )
     &            recip_prs = 1./press0(I,J,bi,bj)
             sig1(I,J) = 0.5*(sigp + sigTmp)*recip_prs
             sig2(I,J) = 0.5*(sigp - sigTmp)*recip_prs
            ENDDO
           ENDDO
           CALL DIAGNOSTICS_FILL(sig1,'SIsig1  ',0,1,2,bi,bj,myThid)
           CALL DIAGNOSTICS_FILL(sig2,'SIsig2  ',0,1,2,bi,bj,myThid)
          ENDDO
         ENDDO
C     endif SEAICEuseEVP
        ENDIF
        DO bj = myByLo(myThid), myByHi(myThid)
         DO bi = myBxLo(myThid), myBxHi(myThid)
          DO j=1,sNy
           DO i=1,sNx
            sigm = e11(I,J,bi,bj) - e22(I,J,bi,bj)
            sigTmp = 
     &           ( e12(I,  J,  bi,bj)**2 + e12(I+1,J,  bi,bj)**2
     &           + e12(I+1,J+1,bi,bj)**2 + e12(I  ,J+1,bi,bj)**2 )
C     shear deformation as sqrt((e11-e22)**2 + 4*e12**2); the 4 pulled into
C     the average
            sig1(I,J) = sqrt(sigm*sigm + sigTmp)
           ENDDO
          ENDDO
          CALL DIAGNOSTICS_FILL(sig1,'SIshear ',0,1,2,bi,bj,myThid)
         ENDDO
        ENDDO
C     endif DIAGNOSTICS_IS_ON(SIsig1/2)
       ENDIF
#endif /* SEAICE_CGRID */
C     abuse press as a temporary field
        IF ( DIAGNOSTICS_IS_ON('SIuheff ',myThid) ) THEN
         DO bj = myByLo(myThid), myByHi(myThid)
          DO bi = myBxLo(myThid), myBxHi(myThid)
           DO j = 1,sNy
            DO i = 1,sNx+1
             press(i,j,bi,bj) =
#ifdef SEAICE_CGRID
     &            UICE(i,j,bi,bj)
#else
C     average B-grid velocities to C-grid points
     &            0.5 _d 0*(UICE(i,j,bi,bj)+UICE(i,j+1,bi,bj))
#endif /* SEAICE_CGRID */
     &            *0.5 _d 0*(HEFF(i,j,bi,bj)+HEFF(i-1,j,bi,bj))
            ENDDO
           ENDDO
          ENDDO
         ENDDO
         CALL DIAGNOSTICS_FILL(press,'SIuheff ',0,1,0,1,1,myThid)
        ENDIF
        IF ( DIAGNOSTICS_IS_ON('SIvheff ',myThid) ) THEN
         DO bj = myByLo(myThid), myByHi(myThid)
          DO bi = myBxLo(myThid), myBxHi(myThid)
           DO j = 1,sNy+1
            DO i = 1,sNx
             press(i,j,bi,bj) =
#ifdef SEAICE_CGRID
     &            VICE(i,j,bi,bj)
#else
C     average B-grid velocities to C-grid points
     &            0.5 _d 0*(VICE(i,j,bi,bj)+VICE(i+1,j,bi,bj))
#endif /* SEAICE_CGRID */
     &            *0.5 _d 0*(HEFF(i,j,bi,bj)+HEFF(i,j-1,bi,bj))
            ENDDO
           ENDDO
          ENDDO
         ENDDO
         CALL DIAGNOSTICS_FILL(press,'SIvheff ',0,1,0,1,1,myThid)
        ENDIF
C     endif useDiagnostics
      ENDIF
#endif /* ALLOW_DIAGNOSTICS */

      RETURN
      END
