C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_init_varia.F,v 1.94 2017/04/04 23:31:27 jmc Exp $
C $Name: $
#include "SEAICE_OPTIONS.h"
#ifdef ALLOW_OBCS
# include "OBCS_OPTIONS.h"
#endif
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 OBCS_UVICE_OLD
# include "OBCS_GRID.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
INTEGER kSurface
_RS mask_uice
INTEGER k
#ifdef ALLOW_SITRACER
INTEGER iTr, jTh
#endif
#ifdef OBCS_UVICE_OLD
INTEGER I_obc, J_obc
#endif /* ALLOW_OBCS */
#ifdef SEAICE_CGRID
_RL recip_tensilDepth
#endif
IF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
kSurface = Nr
ELSE
kSurface = 1
ENDIF
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
HEFFM(i,j,bi,bj) = 0. _d 0
ENDDO
ENDDO
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
HEFFM(i,j,bi,bj)= 1. _d 0
IF (_hFacC(i,j,kSurface,bi,bj).eq.0.)
& HEFFM(i,j,bi,bj)= 0. _d 0
ENDDO
ENDDO
#ifndef SEAICE_CGRID
DO j=1-OLy+1,sNy+OLy
DO i=1-OLx+1,sNx+OLx
UVM(i,j,bi,bj)=0. _d 0
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 _d 0) UVM(i,j,bi,bj)=1. _d 0
ENDDO
ENDDO
#endif /* SEAICE_CGRID */
ENDDO
ENDDO
C coefficients for metric terms
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
#ifdef SEAICE_CGRID
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
k1AtC(I,J,bi,bj) = 0.0 _d 0
k1AtZ(I,J,bi,bj) = 0.0 _d 0
k2AtC(I,J,bi,bj) = 0.0 _d 0
k2AtZ(I,J,bi,bj) = 0.0 _d 0
ENDDO
ENDDO
IF ( usingSphericalPolarGrid .AND. SEAICEuseMetricTerms ) THEN
C This is the only case where tan(phi) is not zero. In this case
C C and U points, and Z and V points have the same phi, so that we
C only need a copy here. Do not use tan(YC) and tan(YG), because
C these
C can be the geographical coordinates and not the correct grid
C coordinates when the grid is rotated (phi/theta/psiEuler .NE. 0)
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
k2AtC(I,J,bi,bj) = - _tanPhiAtU(I,J,bi,bj)*recip_rSphere
k2AtZ(I,J,bi,bj) = - _tanPhiAtV(I,J,bi,bj)*recip_rSphere
ENDDO
ENDDO
ELSEIF ( usingCurvilinearGrid .AND. SEAICEuseMetricTerms ) THEN
C compute metric term coefficients from finite difference
C approximation
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx-1
k1AtC(I,J,bi,bj) = _recip_dyF(I,J,bi,bj)
& * ( _dyG(I+1,J,bi,bj) - _dyG(I,J,bi,bj) )
& * _recip_dxF(I,J,bi,bj)
ENDDO
ENDDO
DO j=1-OLy,sNy+OLy
DO i=1-OLx+1,sNx+OLx
k1AtZ(I,J,bi,bj) = _recip_dyU(I,J,bi,bj)
& * ( _dyC(I,J,bi,bj) - _dyC(I-1,J,bi,bj) )
& * _recip_dxV(I,J,bi,bj)
ENDDO
ENDDO
DO j=1-OLy,sNy+OLy-1
DO i=1-OLx,sNx+OLx
k2AtC(I,J,bi,bj) = _recip_dxF(I,J,bi,bj)
& * ( _dxG(I,J+1,bi,bj) - _dxG(I,J,bi,bj) )
& * _recip_dyF(I,J,bi,bj)
ENDDO
ENDDO
DO j=1-OLy+1,sNy+OLy
DO i=1-OLx,sNx+OLx
k2AtZ(I,J,bi,bj) = _recip_dxV(I,J,bi,bj)
& * ( _dxC(I,J,bi,bj) - _dxC(I,J-1,bi,bj) )
& * _recip_dyU(I,J,bi,bj)
ENDDO
ENDDO
ENDIF
#else /* not SEAICE_CGRID */
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
k1AtC(I,J,bi,bj) = 0.0 _d 0
k1AtU(I,J,bi,bj) = 0.0 _d 0
k1AtV(I,J,bi,bj) = 0.0 _d 0
k2AtC(I,J,bi,bj) = 0.0 _d 0
k2AtU(I,J,bi,bj) = 0.0 _d 0
k2AtV(I,J,bi,bj) = 0.0 _d 0
ENDDO
ENDDO
IF ( usingSphericalPolarGrid .AND. SEAICEuseMetricTerms ) THEN
C This is the only case where tan(phi) is not zero. In this case
C C and U points, and Z and V points have the same phi, so that we
C only need a copy here. Do not use tan(YC) and tan(YG), because
C these
C can be the geographical coordinates and not the correct grid
C coordinates when the grid is rotated (phi/theta/psiEuler .NE. 0)
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
k2AtC(I,J,bi,bj) = - _tanPhiAtU(I,J,bi,bj)*recip_rSphere
k2AtU(I,J,bi,bj) = - _tanPhiAtU(I,J,bi,bj)*recip_rSphere
k2AtV(I,J,bi,bj) = - _tanPhiAtV(I,J,bi,bj)*recip_rSphere
ENDDO
ENDDO
ELSEIF ( usingCurvilinearGrid .AND. SEAICEuseMetricTerms ) THEN
C compute metric term coefficients from finite difference
C approximation
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx-1
k1AtC(I,J,bi,bj) = _recip_dyF(I,J,bi,bj)
& * ( _dyG(I+1,J,bi,bj) - _dyG(I,J,bi,bj) )
& * _recip_dxF(I,J,bi,bj)
ENDDO
ENDDO
DO j=1-OLy,sNy+OLy
DO i=1-OLx+1,sNx+OLx
k1AtU(I,J,bi,bj) = _recip_dyG(I,J,bi,bj)
& * ( _dyF(I,J,bi,bj) - _dyF(I-1,J,bi,bj) )
& * _recip_dxC(I,J,bi,bj)
ENDDO
ENDDO
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx-1
k1AtV(I,J,bi,bj) = _recip_dyC(I,J,bi,bj)
& * ( _dyU(I+1,J,bi,bj) - _dyU(I,J,bi,bj) )
& * _recip_dxG(I,J,bi,bj)
ENDDO
ENDDO
DO j=1-OLy,sNy+OLy-1
DO i=1-OLx,sNx+OLx
k2AtC(I,J,bi,bj) = _recip_dxF(I,J,bi,bj)
& * ( _dxG(I,J+1,bi,bj) - _dxG(I,J,bi,bj) )
& * _recip_dyF(I,J,bi,bj)
ENDDO
ENDDO
DO j=1-OLy,sNy+OLy-1
DO i=1-OLx,sNx+OLx
k2AtU(I,J,bi,bj) = _recip_dxC(I,J,bi,bj)
& * ( _dxV(I,J+1,bi,bj) - _dxV(I,J,bi,bj) )
& * _recip_dyG(I,J,bi,bj)
ENDDO
ENDDO
DO j=1-OLy+1,sNy+OLy
DO i=1-OLx,sNx+OLx
k2AtV(I,J,bi,bj) = _recip_dxG(I,J,bi,bj)
& * ( _dxF(I,J,bi,bj) - _dxF(I,J-1,bi,bj) )
& * _recip_dyC(I,J,bi,bj)
ENDDO
ENDDO
ENDIF
#endif /* not SEAICE_CGRID */
ENDDO
ENDDO
#ifndef SEAICE_CGRID
C-- Choose a proxy level for geostrophic velocity,
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
KGEO(i,j,bi,bj) = 0
ENDDO
ENDDO
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
#ifdef SEAICE_BICE_STRESS
KGEO(i,j,bi,bj) = 1
#else /* SEAICE_BICE_STRESS */
IF (klowc(i,j,bi,bj) .LT. 2) THEN
KGEO(i,j,bi,bj) = 1
ELSE
KGEO(i,j,bi,bj) = 2
DO WHILE ( abs(rC(KGEO(i,j,bi,bj))) .LT. 50.0 _d 0 .AND.
& KGEO(i,j,bi,bj) .LT. (klowc(i,j,bi,bj)-1) )
KGEO(i,j,bi,bj) = KGEO(i,j,bi,bj) + 1
ENDDO
ENDIF
#endif /* SEAICE_BICE_STRESS */
ENDDO
ENDDO
ENDDO
ENDDO
#endif /* SEAICE_CGRID */
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
CToM<<<
#ifdef SEAICE_ITD
DO k=1,nITD
AREAITD(i,j,k,bi,bj) =0. _d 0
HEFFITD(i,j,k,bi,bj) =0. _d 0
ENDDO
#endif
C>>>ToM
UICE(i,j,bi,bj)=0. _d 0
VICE(i,j,bi,bj)=0. _d 0
#ifdef SEAICE_ALLOW_FREEDRIFT
uice_fd(i,j,bi,bj)=0. _d 0
vice_fd(i,j,bi,bj)=0. _d 0
#endif
C
uIceNm1(i,j,bi,bj)=0. _d 0
vIceNm1(i,j,bi,bj)=0. _d 0
ETA (i,j,bi,bj) = 0. _d 0
etaZ(i,j,bi,bj) = 0. _d 0
ZETA(i,j,bi,bj) = 0. _d 0
FORCEX(i,j,bi,bj) = 0. _d 0
FORCEY(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 */
uIceC(i,j,bi,bj) = 0. _d 0
vIceC(i,j,bi,bj) = 0. _d 0
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
#ifdef SEAICE_ALLOW_BOTTOMDRAG
CbotC(i,j,bi,bj) = 0. _d 0
#endif /* SEAICE_ALLOW_BOTTOMDRAG */
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
tensileStrFac(i,j,bi,bj) = 0. _d 0
CToM<<<
#ifdef SEAICE_ITD
DO k=1,nITD
HSNOWITD(i,j,k,bi,bj)=0. _d 0
ENDDO
#endif
C>>>ToM
#ifdef SEAICE_VARIABLE_SALINITY
HSALT(i,j,bi,bj) = 0. _d 0
#endif
#ifdef ALLOW_SITRACER
DO iTr = 1, SItrMaxNum
SItracer(i,j,bi,bj,iTr) = 0. _d 0
SItrBucket(i,j,bi,bj,iTr) = 0. _d 0
c "ice concentration" tracer that should remain .EQ.1.
if (SItrName(iTr).EQ.'one') SItracer(i,j,bi,bj,iTr)=1. _d 0
ENDDO
DO jTh = 1, 5
SItrHEFF (i,j,bi,bj,jTh) = 0. _d 0
ENDDO
DO jTh = 1, 3
SItrAREA (i,j,bi,bj,jTh) = 0. _d 0
ENDDO
#endif
DO k=1,nITD
TICES(i,j,k,bi,bj)=0. _d 0
ENDDO
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
icevolMeanCell(i,j,bi,bj) = 0. _d 0
#endif
saltWtrIce(i,j,bi,bj) = 0. _d 0
frWtrIce(i,j,bi,bj) = 0. _d 0
#if (defined (ALLOW_MEAN_SFLUX_COST_CONTRIBUTION) defined (ALLOW_SSH_GLOBMEAN_COST_CONTRIBUTION))
frWtrAtm(i,j,bi,bj) = 0. _d 0
AREAforAtmFW(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
#ifdef SEAICE_CGRID
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
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
ENDDO
ENDDO
ENDDO
ENDDO
#endif /* SEAICE_CGRID */
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
#ifdef OBCS_UVICE_OLD
#ifdef SEAICE_CGRID
IF (useOBCS) THEN
C-- If OBCS is turned on, close southern and western boundaries
DO i=1-OLx,sNx+OLx
C Southern boundary
J_obc = OB_Js(i,bi,bj)
IF ( J_obc.NE.OB_indexNone ) THEN
seaiceMaskU(i,J_obc,bi,bj)= 0.0 _d 0
seaiceMaskV(i,J_obc,bi,bj)= 0.0 _d 0
ENDIF
ENDDO
DO j=1-OLy,sNy+OLy
C Western boundary
I_obc = OB_Iw(j,bi,bj)
IF ( I_obc.NE.OB_indexNone ) THEN
seaiceMaskU(I_obc,j,bi,bj)= 0.0 _d 0
seaiceMaskV(I_obc,j,bi,bj)= 0.0 _d 0
ENDIF
ENDDO
ENDIF
#endif /* SEAICE_CGRID */
#endif /* OBCS_UVICE_OLD */
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
DO k=1,nITD
TICES(i,j,k,bi,bj)=273.0 _d 0
ENDDO
#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_RS(UVM, myThid)
#endif
C-- Now lets look at all these beasts
IF ( plotLevel.GE.debLevC ) THEN
CALL PLOT_FIELD_XYRL( HEFFM , 'Current HEFFM ' ,
& nIter0, myThid )
#ifdef SEAICE_CGRID
CALL PLOT_FIELD_XYRL( seaiceMaskU, 'Current seaiceMaskU',
& nIter0, myThid )
CALL PLOT_FIELD_XYRL( seaiceMaskV, 'Current seaiceMaskV',
& nIter0, myThid )
#else
CALL PLOT_FIELD_XYRS( UVM , 'Current UVM ' ,
& nIter0, 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
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 velocity from file (if available)
IF ( uIceFile .NE. ' ' )
& CALL READ_FLD_XY_RL( uIceFile, ' ', uIce, 0, myThid )
IF ( vIceFile .NE. ' ' )
& CALL READ_FLD_XY_RL( vIceFile, ' ', vIce, 0, myThid )
IF ( uIceFile .NE. ' ' .OR. vIceFile .NE. ' ' ) THEN
#ifdef SEAICE_CGRID
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
uIce(i,j,bi,bj) = uIce(i,j,bi,bj)*seaiceMaskU(i,j,bi,bj)
vIce(i,j,bi,bj) = vIce(i,j,bi,bj)*seaiceMaskV(i,j,bi,bj)
ENDDO
ENDDO
ENDDO
ENDDO
#endif /* SEAICE_CGRID */
CALL EXCH_UV_XY_RL( uIce, vIce, .TRUE., myThid )
ENDIF
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_ITD
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
AREAITD(I,J,1,bi,bj) = AREA(I,J,bi,bj)
HEFFITD(I,J,1,bi,bj) = HEFF(I,J,bi,bj)
HSNOWITD(I,J,1,bi,bj) = HSNOW(I,J,bi,bj)
opnWtrFrac(I,J,bi,bj) = 1. _d 0 - AREA(I,J,bi,bj)
fw2ObyRidge(I,J,bi,bj) = 0. _d 0
ENDDO
ENDDO
CALL SEAICE_ITD_REDIST(bi, bj, baseTime, nIter0, myThid)
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*SEAICE_saltFrac
cif & ICE2WATR*rhoConstFresh*SEAICE_saltFrac
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 ALLOW_SITRACER
C-- Read initial sea ice age from file if available.
DO iTr = 1, SItrMaxNum
IF ( SItrFile(iTr) .NE. ' ' ) THEN
CALL READ_FLD_XY_RL( siTrFile(iTr), ' ',
& SItracer(1-OLx,1-OLy,1,1,iTr), 0, myThid )
_EXCH_XY_RL(SItracer(1-OLx,1-OLy,1,1,iTr),myThid)
ENDIF
ENDDO
#endif /* ALLOW_SITRACER */
ENDIF
#if (defined (ALLOW_MEAN_SFLUX_COST_CONTRIBUTION) defined (ALLOW_SSH_GLOBMEAN_COST_CONTRIBUTION))
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
AREAforAtmFW(i,j,bi,bj) = AREA(i,j,bi,bj)
ENDDO
ENDDO
ENDDO
ENDDO
#endif
#ifdef ALLOW_OBCS
C-- In case we use scheme with a large stencil that extends into overlap:
C no longer needed with the right masking in advection & diffusion S/R.
c IF ( useOBCS ) THEN
c DO bj=myByLo(myThid),myByHi(myThid)
c DO bi=myBxLo(myThid),myBxHi(myThid)
c CALL OBCS_COPY_TRACER( HEFF(1-OLx,1-OLy,bi,bj),
c I 1, bi, bj, myThid )
c CALL OBCS_COPY_TRACER( AREA(1-OLx,1-OLy,bi,bj),
c I 1, bi, bj, myThid )
c CALL OBCS_COPY_TRACER( HSNOW(1-OLx,1-OLy,bi,bj),
c I 1, bi, bj, myThid )
#ifdef SEAICE_VARIABLE_SALINITY
c CALL OBCS_COPY_TRACER( HSALT(1-OLx,1-OLy,bi,bj),
c I 1, bi, bj, myThid )
#endif
c ENDDO
c ENDDO
c ENDIF
#endif /* ALLOW_OBCS */
#ifdef SEAICE_ALLOW_JFNK
C Computing this metric cannot be done in S/R SEAICE_INIT_FIXED
C where it belongs, because globalArea is only defined later after
C S/R PACKAGES_INIT_FIXED, so we move this computation here.
CALL SEAICE_MAP_RS2VEC( nVec, rAw, rAs,
& scalarProductMetric, .TRUE., myThid )
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
DO i=1,nVec
scalarProductMetric(i,1,bi,bj) =
& scalarProductMetric(i,1,bi,bj)/globalArea
ENDDO
ENDDO
ENDDO
#endif /* SEAICE_ALLOW_JFNK */
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,bi,bj)*(1.0 _d 11)
ETA(i,j,bi,bj) = ZETA(i,j,bi,bj)/SEAICE_eccen**2
PRESS0(i,j,bi,bj) = SEAICE_strength*HEFF(i,j,bi,bj)
& *EXP(-SEAICE_cStar*(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
#ifdef SEAICE_CGRID
C compute tensile strength factor k: tensileStrength = k*PRESS
C can be done in initialisation phase as long as it depends only
C on depth
IF ( SEAICE_tensilFac .NE. 0. _d 0 ) THEN
recip_tensilDepth = 0. _d 0
IF ( SEAICE_tensilDepth .GT. 0. _d 0 )
& recip_tensilDepth = 1. _d 0 / SEAICE_tensilDepth
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
tensileStrFac(i,j,bi,bj) = SEAICE_tensilFac*HEFFM(i,j,bi,bj)
& *exp(-ABS(R_low(i,j,bi,bj))*recip_tensilDepth)
ENDDO
ENDDO
ENDDO
ENDDO
ENDIF
#endif /* SEAICE_CGRID */
RETURN
END