C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_init_fixed.F,v 1.30 2014/08/23 20:22:16 torge Exp $
C $Name:  $

#include "SEAICE_OPTIONS.h"

CStartOfInterface
      SUBROUTINE SEAICE_INIT_FIXED( myThid )
C     *==========================================================*
C     | SUBROUTINE SEAICE_INIT_FIXED
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 "FFIELDS.h"
#include "SEAICE_SIZE.h"
#include "SEAICE_PARAMS.h"
#include "SEAICE.h"
#include "SEAICE_TRACER.h"

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

C     === Local variables ===
#ifdef SEAICE_ITD
C     msgBuf      :: Informational/error message buffer
      CHARACTER*(MAX_LEN_MBUF) msgBuf
      CHARACTER*15 HlimitMsgFormat
C     k - loop counter for ITD categories
      INTEGER k
      _RL tmpVar
#endif
C     i,j, bi,bj  :: Loop counters
c     INTEGER i, j, bi, bj
      INTEGER kSurface
#ifdef ALLOW_SITRACER
      INTEGER iTracer
#endif
#ifdef SHORTWAVE_HEATING
cif   Helper variable for determining the fraction of sw radiation
cif   penetrating the model shallowest layer
      _RL dummyTime
      _RL swfracba(2)
      _RL tmpFac
#endif /* SHORTWAVE_HEATING */

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

C     Initialize MNC variable information for SEAICE
      IF ( useMNC .AND.
     &    (seaice_tave_mnc.OR.seaice_dump_mnc.OR.SEAICE_mon_mnc)
     &   ) THEN
        CALL SEAICE_MNC_INIT( myThid )
      ENDIF

C     Only Master Thread updates parameter in commom block:
      _BEGIN_MASTER(myThid)

C     restart parameter
      SEAICEmomStartBDF = 0
      IF ( SEAICEuseBDF2 ) SEAICEmomStartBDF = nIter0

#ifdef SHORTWAVE_HEATING
       tmpFac    = -1.0
       dummyTime = 1.0
       swfracba(1) = ABS(rF(1))
       swfracba(2) = ABS(rF(2))
       CALL SWFRAC(
     I       2, tmpFac,
     U       swfracba,
     I       dummyTime, 0, myThid )
       SWFracB = swfracba(2)
#else /* SHORTWAVE_HEATING */
       SWFracB = 0. _d 0
#endif /* SHORTWAVE_HEATING */

C--   Set mcPheePiston coeff (if still unset)
      IF ( SEAICE_mcPheePiston.EQ.UNSET_RL ) THEN
        IF ( SEAICE_availHeatFrac.NE.UNSET_RL ) THEN
          SEAICE_mcPheePiston = SEAICE_availHeatFrac
     &                        * drF(kSurface)/SEAICE_deltaTtherm
        ELSE
          SEAICE_mcPheePiston = MCPHEE_TAPER_FAC
     &                        * STANTON_NUMBER * USTAR_BASE
          SEAICE_mcPheePiston = MIN( SEAICE_mcPheePiston,
     &                          drF(kSurface)/SEAICE_deltaTtherm )
        ENDIF
      ENDIF

C--   SItracer specifications for basic tracers
#ifdef ALLOW_SITRACER
      DO iTracer = 1, SItrNumInUse
C     "ice concentration" tracer that should remain .EQ.1.
       IF (SItrName(iTracer).EQ.'one') THEN
         SItrFromOcean0(iTracer)    =ONE
         SItrFromFlood0(iTracer)    =ONE
         SItrExpand0(iTracer)       =ONE
         SItrFromOceanFrac(iTracer) =ZERO
         SItrFromFloodFrac(iTracer) =ZERO
       ENDIF
C     age tracer: no age in ocean, or effect from ice cover changes
       IF (SItrName(iTracer).EQ.'age') THEN
         SItrFromOcean0(iTracer)    =ZERO
         SItrFromFlood0(iTracer)    =ZERO
         SItrExpand0(iTracer)       =ZERO
         SItrFromOceanFrac(iTracer) =ZERO
         SItrFromFloodFrac(iTracer) =ZERO
       ENDIf
C     salinity tracer:
       IF (SItrName(iTracer).EQ.'salinity') THEN
         SItrMate(iTracer)          ='HEFF'
         SItrExpand0(iTracer)       =ZERO
         IF ( SEAICE_salinityTracer ) THEN
           SEAICE_salt0    = ZERO
           SEAICE_saltFrac = ZERO
         ENDIF
       ENDIF
C     simple, made up, ice surface roughness index prototype
       IF (SItrName(iTracer).EQ.'ridge') THEN
         SItrMate(iTracer)          ='AREA'
         SItrFromOcean0(iTracer)    =ZERO
         SItrFromFlood0(iTracer)    =ZERO
         SItrExpand0(iTracer)       =ZERO
         SItrFromOceanFrac(iTracer) =ZERO
         SItrFromFloodFrac(iTracer) =ZERO
       ENDIF
#ifdef SEAICE_GREASE
C     grease ice tracer:
c     (Smedrud and Martin, 2014, Ann. Glac.)
       IF (SItrName(iTracer).EQ.'grease') THEN
         SItrMate(iTracer)          ='HEFF'
         SItrFromOcean0(iTracer)    =ZERO
         SItrFromFlood0(iTracer)    =ZERO
         SItrExpand0(iTracer)       =ZERO
         SItrFromOceanFrac(iTracer) =ZERO
         SItrFromFloodFrac(iTracer) =ZERO
       ENDIF
#endif /* SEAICE_GREASE */
      ENDDO
#endif /* ALLOW_SITRACER */

#ifdef SEAICE_ITD
C     use Equ. 22 of Lipscomb et al. (2001, JGR) to generate ice
C     thickness category limits:
C     - dependends on given number of categories nITD
C     - choose between original parameters of Lipscomb et al. (2001):
C       c1=3.0/N, c2=15*c1, c3=3.0
C       or emphasize thin end of ITD (in order to enhance ice growth):
C       c1=1.5/N, c2=42*c1, c3=3.3
C       -> HINT: set parameters c1, c2 and c3 in seaice_readparms.F
      Hlimit(0) = 0. _d 0
      IF ( nITD.GT.1 ) THEN
       tmpVar = nITD
       tmpVar = oneRL / tmpVar
       Hlimit_c1 = Hlimit_c1*tmpVar
       Hlimit_c2 = Hlimit_c2*Hlimit_c1
       DO k=1,nITD-1
        Hlimit(k) = Hlimit(k-1)
     &            + Hlimit_c1
     &            + Hlimit_c2
     &  *( oneRL + TANH( Hlimit_c3 *( FLOAT(k-1)*tmpVar - oneRL ) ) )
       ENDDO
      ENDIF
C     thickest category is unlimited
      Hlimit(nITD) = 999.9 _d 0

      WRITE(msgBuf,'(A,I2,A)')
     &     ' SEAICE_INIT_FIXED: ', nITD,
     &     ' sea ice thickness categories'
      CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                    SQUEEZE_RIGHT, myThid )
      WRITE(HlimitMsgFormat,'(A,I2,A)') '(A,',nITD,'F6.2,F6.1)'
      WRITE(msgBuf,HlimitMsgFormat)
     &     ' SEAICE_INIT_FIXED: Hlimit = ', (Hlimit(k),k=0,nITD)
      CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                    SQUEEZE_RIGHT, myThid )
#endif /* SEAICE_ITD */

C     Only Master Thread updates parameter in common block:
      _END_MASTER(myThid)

#ifdef SEAICE_ALLOW_JFNK
C     Only Master Thread updates parameter in common block:
      _BEGIN_MASTER(myThid)
C     initialise some diagnostic counters for the JFNK solver
      totalNewtonIters   = 0
      totalNewtonFails   = 0
      totalKrylovIters   = 0
      totalKrylovFails   = 0
      totalJFNKtimeSteps = 0
      _END_MASTER(myThid)
C     this cannot be done here, because globalArea is only defined
C     after S/R PACKAGES_INIT_FIXED, so we move it to S/R SEAICE_INIT_VARIA
CML      CALL SEAICE_MAP2VEC( nVec, rAw, rAs,
CML     &     scalarProductMetric, .TRUE., myThid )
CML      DO bj=myByLo(myThid),myByHi(myThid)
CML       DO bi=myBxLo(myThid),myBxHi(myThid)
CML        DO i=1,nVec
CML         scalarProductMetric(i,1,bi,bj) =
CML     &        scalarProductMetric(i,1,bi,bj)/globalArea
CML        ENDDO
CML       ENDDO
CML      ENDDO
#endif /* SEAICE_ALLOW_JFNK */

C--   all threads wait for master to finish initialisation of shared params
      _BARRIER

#ifdef ALLOW_DIAGNOSTICS
      IF ( useDiagnostics ) THEN
        CALL SEAICE_DIAGNOSTICS_INIT( myThid )
      ENDIF
#endif

C--   Summarise pkg/seaice configuration
      CALL SEAICE_SUMMARY( myThid )

      RETURN
      END