Back to home page

MITgcm

 
 

    


File indexing completed on 2025-03-03 06:11:47 UTC

view on githubraw file Latest commit b7b61e61 on 2025-03-02 15:55:22 UTC
809c36b928 Patr*0001 #include "SEAICE_OPTIONS.h"
fc004572cd Jean*0002 #ifdef ALLOW_OBCS
                0003 # include "OBCS_OPTIONS.h"
                0004 #endif
809c36b928 Patr*0005 
                0006 CStartOfInterface
7e4f7741f6 Patr*0007       SUBROUTINE SEAICE_INIT_VARIA( myThid )
e5b8ebbabe Jean*0008 C     *==========================================================*
7e4f7741f6 Patr*0009 C     | SUBROUTINE SEAICE_INIT_VARIA                             |
809c36b928 Patr*0010 C     | o Initialization of sea ice model.                       |
e5b8ebbabe Jean*0011 C     *==========================================================*
                0012 C     *==========================================================*
809c36b928 Patr*0013       IMPLICIT NONE
f4cb8ca7db Jean*0014 
809c36b928 Patr*0015 C     === Global variables ===
                0016 #include "SIZE.h"
                0017 #include "EEPARAMS.h"
                0018 #include "PARAMS.h"
                0019 #include "GRID.h"
36c84b38c7 Dimi*0020 #include "DYNVARS.h"
f3ce416a61 Jean*0021 #include "FFIELDS.h"
ccaa3c61f4 Patr*0022 #include "SEAICE_SIZE.h"
8c50aa8796 Jean*0023 #include "SEAICE_PARAMS.h"
                0024 #include "SEAICE.h"
ccaa3c61f4 Patr*0025 #include "SEAICE_TRACER.h"
8c50aa8796 Jean*0026 #include "SEAICE_TAVE.h"
fc004572cd Jean*0027 #ifdef OBCS_UVICE_OLD
a9eb030de2 Jean*0028 # include "OBCS_GRID.h"
ae1fb66b64 Dimi*0029 #endif
809c36b928 Patr*0030 
                0031 C     === Routine arguments ===
425e8efc36 Jean*0032 C     myThid :: Thread no. that called this routine.
809c36b928 Patr*0033       INTEGER myThid
                0034 CEndOfInterface
f4cb8ca7db Jean*0035 
809c36b928 Patr*0036 C     === Local variables ===
425e8efc36 Jean*0037 C     i,j,k,bi,bj :: Loop counters
809c36b928 Patr*0038 
edfdf5fa1d Jean*0039       INTEGER i, j, bi, bj
0320e25227 Mart*0040       INTEGER kSrf
8c50aa8796 Jean*0041       INTEGER k
f50f58ec54 Gael*0042 #ifdef ALLOW_SITRACER
                0043       INTEGER iTr, jTh
                0044 #endif
fc004572cd Jean*0045 #ifdef OBCS_UVICE_OLD
dc8f395c29 Dimi*0046       INTEGER I_obc, J_obc
                0047 #endif /* ALLOW_OBCS */
0dcda34642 Jean*0048 #ifdef SEAICE_CGRID
45315406aa Mart*0049       _RS mask_uice
2f5e8addfd Mart*0050       _RL recip_tensilDepth
0dcda34642 Jean*0051 #endif
486cba1d02 Mart*0052 
0320e25227 Mart*0053       IF ( usingPCoords ) THEN
                0054        kSrf = Nr
f4cb8ca7db Jean*0055       ELSE
0320e25227 Mart*0056        kSrf = 1
f4cb8ca7db Jean*0057       ENDIF
a6d2313c04 Dimi*0058 
f4cb8ca7db Jean*0059 C--   Initialise all variables in common blocks:
7109a141b2 Patr*0060       DO bj=myByLo(myThid),myByHi(myThid)
                0061        DO bi=myBxLo(myThid),myBxHi(myThid)
04016f2c47 Mart*0062         DO j=1-OLy,sNy+OLy
                0063          DO i=1-OLx,sNx+OLx
772590b63c Mart*0064           HEFF(i,j,bi,bj)=0. _d 0
                0065           AREA(i,j,bi,bj)=0. _d 0
45315406aa Mart*0066           HSNOW(i,j,bi,bj)   = 0. _d 0
86b84a92fc Patr*0067 #ifdef SEAICE_ITD
                0068           DO k=1,nITD
                0069            AREAITD(i,j,k,bi,bj) =0. _d 0
                0070            HEFFITD(i,j,k,bi,bj) =0. _d 0
45315406aa Mart*0071            HSNOWITD(i,j,k,bi,bj)=0. _d 0
86b84a92fc Patr*0072           ENDDO
                0073 #endif
772590b63c Mart*0074           UICE(i,j,bi,bj)=0. _d 0
                0075           VICE(i,j,bi,bj)=0. _d 0
                0076 C
5fe78992ba Mart*0077           uIceNm1(i,j,bi,bj) = 0. _d 0
                0078           vIceNm1(i,j,bi,bj) = 0. _d 0
45315406aa Mart*0079           DWATN  (i,j,bi,bj) = 0. _d 0
                0080 #ifdef SEAICE_CGRID
                0081           stressDivergenceX(i,j,bi,bj) = 0. _d 0
                0082           stressDivergenceY(i,j,bi,bj) = 0. _d 0
                0083 # ifdef SEAICE_ALLOW_EVP
                0084           seaice_sigma1 (i,j,bi,bj) = 0. _d 0
                0085           seaice_sigma2 (i,j,bi,bj) = 0. _d 0
                0086           seaice_sigma12(i,j,bi,bj) = 0. _d 0
                0087 # endif /* SEAICE_ALLOW_EVP */
                0088 #endif
                0089 #if ( defined SEAICE_CGRID || defined SEAICE_BGRID_DYNAMICS )
5fe78992ba Mart*0090           e11    (i,j,bi,bj) = 0. _d 0
                0091           e22    (i,j,bi,bj) = 0. _d 0
                0092           e12    (i,j,bi,bj) = 0. _d 0
                0093           deltaC (i,j,bi,bj) = 0. _d 0
45315406aa Mart*0094           PRESS  (i,j,bi,bj) = 0. _d 0
5fe78992ba Mart*0095           ETA    (i,j,bi,bj) = 0. _d 0
                0096           etaZ   (i,j,bi,bj) = 0. _d 0
                0097           ZETA   (i,j,bi,bj) = 0. _d 0
                0098           FORCEX (i,j,bi,bj) = 0. _d 0
                0099           FORCEY (i,j,bi,bj) = 0. _d 0
45315406aa Mart*0100           tensileStrFac(i,j,bi,bj) = 0. _d 0
                0101           PRESS0(i,j,bi,bj)  = 0. _d 0
                0102           FORCEX0(i,j,bi,bj) = 0. _d 0
                0103           FORCEY0(i,j,bi,bj) = 0. _d 0
                0104           SEAICE_zMax(i,j,bi,bj) = 0. _d 0
                0105           SEAICE_zMin(i,j,bi,bj) = 0. _d 0
                0106 #endif
f4cb8ca7db Jean*0107 #ifdef SEAICE_CGRID
                0108           seaiceMassC(i,j,bi,bj)=0. _d 0
                0109           seaiceMassU(i,j,bi,bj)=0. _d 0
                0110           seaiceMassV(i,j,bi,bj)=0. _d 0
45315406aa Mart*0111 # ifdef SEAICE_ALLOW_FREEDRIFT
                0112           uice_fd(i,j,bi,bj) = 0. _d 0
                0113           vice_fd(i,j,bi,bj) = 0. _d 0
                0114 # endif
                0115 # ifdef SEAICE_ALLOW_BOTTOMDRAG
                0116           CbotC(i,j,bi,bj)   = 0. _d 0
                0117 # endif /* SEAICE_ALLOW_BOTTOMDRAG */
5bb179ddc2 Mart*0118 # ifdef SEAICE_ALLOW_SIDEDRAG
                0119           sideDragU(i,j,bi,bj) = 0. _d 0
                0120           sideDragV(i,j,bi,bj) = 0. _d 0
                0121 # endif /* SEAICE_ALLOW_SIDEDRAG */
45315406aa Mart*0122 #endif /* SEAICE_CGRID */
                0123 #ifdef SEAICE_BGRID_DYNAMICS
                0124           uIceB(i,j,bi,bj)   = 0. _d 0
                0125           vIceB(i,j,bi,bj)   = 0. _d 0
5fe78992ba Mart*0126           AMASS(i,j,bi,bj)   = 0. _d 0
                0127           DAIRN(i,j,bi,bj)   = 0. _d 0
                0128           WINDX(i,j,bi,bj)   = 0. _d 0
                0129           WINDY(i,j,bi,bj)   = 0. _d 0
                0130           GWATX(i,j,bi,bj)   = 0. _d 0
                0131           GWATY(i,j,bi,bj)   = 0. _d 0
45315406aa Mart*0132 #endif /* SEAICE_BGRID_DYNAMICS */
a98c4b8072 Ian *0133 #ifdef SEAICE_VARIABLE_SALINITY
f4cb8ca7db Jean*0134           HSALT(i,j,bi,bj)  = 0. _d 0
                0135 #endif
f50f58ec54 Gael*0136 #ifdef ALLOW_SITRACER
                0137           DO iTr = 1, SItrMaxNum
                0138            SItracer(i,j,bi,bj,iTr) = 0. _d 0
                0139            SItrBucket(i,j,bi,bj,iTr) = 0. _d 0
                0140 c "ice concentration" tracer that should remain .EQ.1.
                0141            if (SItrName(iTr).EQ.'one') SItracer(i,j,bi,bj,iTr)=1. _d 0
                0142           ENDDO
                0143           DO jTh = 1, 5
                0144            SItrHEFF (i,j,bi,bj,jTh) = 0. _d 0
                0145           ENDDO
bb24b8a3e6 Gael*0146           DO jTh = 1, 3
                0147            SItrAREA (i,j,bi,bj,jTh) = 0. _d 0
                0148           ENDDO
f50f58ec54 Gael*0149 #endif
f913c5a485 Mart*0150           DO k=1,nITD
f5282c5b03 Gael*0151             TICES(i,j,k,bi,bj)=0. _d 0
                0152           ENDDO
1cfef3b6ad Patr*0153           saltWtrIce(i,j,bi,bj) = 0. _d 0
                0154           frWtrIce(i,j,bi,bj)   = 0. _d 0
7109a141b2 Patr*0155          ENDDO
                0156         ENDDO
                0157        ENDDO
                0158       ENDDO
                0159 
c69ccc91fb antn*0160 #ifdef ALLOW_AUTODIFF
                0161 C- Note: To simplify dependency & avoid recomputations, when compiling
                0162 C        pkg/autodiff, we always call SEAICE_INIT_VARIA to initialise control
                0163 C        variables (as done above) without condition on useSEAICE.
                0164 C        Therefore, in this case, the "If useSEAICE" is added back here:
                0165       IF ( useSEAICE ) THEN
                0166 #endif
                0167 
8c50aa8796 Jean*0168 #ifdef ALLOW_TIMEAVE
                0169 C     Initialize averages to zero
                0170       DO bj = myByLo(myThid), myByHi(myThid)
                0171        DO bi = myBxLo(myThid), myBxHi(myThid)
                0172         CALL TIMEAVE_RESET( FUtave   , 1, bi, bj, myThid )
                0173         CALL TIMEAVE_RESET( FVtave   , 1, bi, bj, myThid )
                0174         CALL TIMEAVE_RESET( EmPmRtave, 1, bi, bj, myThid )
                0175         CALL TIMEAVE_RESET( QNETtave , 1, bi, bj, myThid )
                0176         CALL TIMEAVE_RESET( QSWtave  , 1, bi, bj, myThid )
                0177         CALL TIMEAVE_RESET( UICEtave , 1, bi, bj, myThid )
                0178         CALL TIMEAVE_RESET( VICEtave , 1, bi, bj, myThid )
                0179         CALL TIMEAVE_RESET( HEFFtave , 1, bi, bj, myThid )
                0180         CALL TIMEAVE_RESET( AREAtave , 1, bi, bj, myThid )
                0181         SEAICE_timeAve(bi,bj) = ZERO
                0182        ENDDO
                0183       ENDDO
                0184 #endif /* ALLOW_TIMEAVE */
                0185 
84cb676f82 Mart*0186 C--   Initialize (variable) grid info. As long as we allow masking of
                0187 C--   velocities outside of ice covered areas (in seaice_dynsolver)
b7b61e618a Mart*0188 C--   we need to re-initialize seaiceMaskU/V here for TAF
809c36b928 Patr*0189       DO bj=myByLo(myThid),myByHi(myThid)
                0190        DO bi=myBxLo(myThid),myBxHi(myThid)
45315406aa Mart*0191 #ifdef SEAICE_CGRID
f4cb8ca7db Jean*0192         DO j=1-OLy+1,sNy+OLy
                0193          DO i=1-OLx+1,sNx+OLx
                0194           seaiceMaskU(i,j,bi,bj)=   0.0 _d 0
                0195           seaiceMaskV(i,j,bi,bj)=   0.0 _d 0
                0196           mask_uice=HEFFM(i,j,bi,bj)+HEFFM(i-1,j  ,bi,bj)
a808b72b4b Patr*0197           IF(mask_uice.GT.1.5 _d 0) seaiceMaskU(i,j,bi,bj)=1.0 _d 0
f4cb8ca7db Jean*0198           mask_uice=HEFFM(i,j,bi,bj)+HEFFM(i  ,j-1,bi,bj)
a808b72b4b Patr*0199           IF(mask_uice.GT.1.5 _d 0) seaiceMaskV(i,j,bi,bj)=1.0 _d 0
ab37184352 Mart*0200          ENDDO
                0201         ENDDO
45315406aa Mart*0202 # ifdef OBCS_UVICE_OLD
eb0337bb06 Dimi*0203         IF (useOBCS) THEN
b48f786726 Dimi*0204 C--   If OBCS is turned on, close southern and western boundaries
48474b8e08 Jean*0205          DO i=1-OLx,sNx+OLx
a6d2313c04 Dimi*0206 C Southern boundary
dc8f395c29 Dimi*0207           J_obc = OB_Js(i,bi,bj)
d6145066cc Jean*0208           IF ( J_obc.NE.OB_indexNone ) THEN
dc8f395c29 Dimi*0209            seaiceMaskU(i,J_obc,bi,bj)=   0.0 _d 0
                0210            seaiceMaskV(i,J_obc,bi,bj)=   0.0 _d 0
eb0337bb06 Dimi*0211           ENDIF
                0212          ENDDO
48474b8e08 Jean*0213          DO j=1-OLy,sNy+OLy
a6d2313c04 Dimi*0214 C Western boundary
d6145066cc Jean*0215           I_obc = OB_Iw(j,bi,bj)
                0216           IF ( I_obc.NE.OB_indexNone ) THEN
dc8f395c29 Dimi*0217            seaiceMaskU(I_obc,j,bi,bj)=   0.0 _d 0
                0218            seaiceMaskV(I_obc,j,bi,bj)=   0.0 _d 0
b48f786726 Dimi*0219           ENDIF
809c36b928 Patr*0220          ENDDO
b48f786726 Dimi*0221         ENDIF
45315406aa Mart*0222 # endif /* OBCS_UVICE_OLD */
a6d2313c04 Dimi*0223 #endif /* SEAICE_CGRID */
69c991a1ec Dimi*0224 
809c36b928 Patr*0225         DO j=1-OLy,sNy+OLy
                0226          DO i=1-OLx,sNx+OLx
f913c5a485 Mart*0227           DO k=1,nITD
f4cb8ca7db Jean*0228            TICES(i,j,k,bi,bj)=273.0 _d 0
09510da3bb Dimi*0229           ENDDO
45315406aa Mart*0230 #ifdef SEAICE_CGRID
f4cb8ca7db Jean*0231           seaiceMassC(i,j,bi,bj)=1000.0 _d 0
                0232           seaiceMassU(i,j,bi,bj)=1000.0 _d 0
                0233           seaiceMassV(i,j,bi,bj)=1000.0 _d 0
53092bcb42 Mart*0234 #endif
45315406aa Mart*0235 #ifdef SEAICE_BGRID_DYNAMICS
                0236           AMASS      (i,j,bi,bj)=1000.0 _d 0
                0237 #endif
809c36b928 Patr*0238          ENDDO
                0239         ENDDO
d5bc48d0a9 Dimi*0240 
809c36b928 Patr*0241        ENDDO
                0242       ENDDO
                0243 
                0244 C--   Update overlap regions
3ed8349f04 Mart*0245 #ifdef SEAICE_CGRID
b4f21203f7 Jean*0246       CALL EXCH_UV_XY_RL(seaiceMaskU,seaiceMaskV,.FALSE.,myThid)
45315406aa Mart*0247 #endif
                0248 #ifdef SEAICE_BGRID_DYNAMICS
48474b8e08 Jean*0249       _EXCH_XY_RS(UVM, myThid)
3ed8349f04 Mart*0250 #endif
809c36b928 Patr*0251 
a9a8db28c6 Alis*0252 C--   Now lets look at all these beasts
a24915ab1a Jean*0253       IF ( plotLevel.GE.debLevC ) THEN
3e0af136db Dimi*0254          CALL PLOT_FIELD_XYRL( HEFFM   , 'Current HEFFM   ' ,
23142459d0 Jean*0255      &        nIter0, myThid )
3ed8349f04 Mart*0256 #ifdef SEAICE_CGRID
                0257          CALL PLOT_FIELD_XYRL( seaiceMaskU, 'Current seaiceMaskU',
23142459d0 Jean*0258      &        nIter0, myThid )
3ed8349f04 Mart*0259          CALL PLOT_FIELD_XYRL( seaiceMaskV, 'Current seaiceMaskV',
23142459d0 Jean*0260      &        nIter0, myThid )
45315406aa Mart*0261 #endif
                0262 #ifdef SEAICE_BGRID_DYNAMICS
48474b8e08 Jean*0263          CALL PLOT_FIELD_XYRS( UVM     , 'Current UVM     ' ,
23142459d0 Jean*0264      &        nIter0, myThid )
3ed8349f04 Mart*0265 #endif
3e0af136db Dimi*0266       ENDIF
09510da3bb Dimi*0267 
809c36b928 Patr*0268 C--   Set model variables to initial/restart conditions
2a682dc462 Mart*0269       IF ( .NOT. ( startTime .EQ. baseTime .AND.  nIter0 .EQ. 0
                0270      &     .AND. pickupSuff .EQ. ' ') ) THEN
36f6faaf99 Dimi*0271 
6060ec2938 Dimi*0272          CALL SEAICE_READ_PICKUP ( myThid )
36f6faaf99 Dimi*0273 
6060ec2938 Dimi*0274       ELSE
36f6faaf99 Dimi*0275 
e7c33124ce Dimi*0276        DO bj=myByLo(myThid),myByHi(myThid)
                0277         DO bi=myBxLo(myThid),myBxHi(myThid)
                0278          DO j=1-OLy,sNy+OLy
                0279           DO i=1-OLx,sNx+OLx
772590b63c Mart*0280            HEFF(i,j,bi,bj)=SEAICE_initialHEFF*HEFFM(i,j,bi,bj)
                0281            UICE(i,j,bi,bj)=ZERO
                0282            VICE(i,j,bi,bj)=ZERO
e7c33124ce Dimi*0283           ENDDO
809c36b928 Patr*0284          ENDDO
e7c33124ce Dimi*0285         ENDDO
                0286        ENDDO
6060ec2938 Dimi*0287 
425e8efc36 Jean*0288 C--   Read initial sea-ice velocity from file (if available)
                0289        IF ( uIceFile .NE. ' ' )
                0290      &  CALL READ_FLD_XY_RL( uIceFile, ' ', uIce, 0, myThid )
                0291        IF ( vIceFile .NE. ' ' )
                0292      &  CALL READ_FLD_XY_RL( vIceFile, ' ', vIce, 0, myThid )
                0293        IF ( uIceFile .NE. ' ' .OR. vIceFile .NE. ' ' ) THEN
                0294 #ifdef SEAICE_CGRID
                0295         DO bj=myByLo(myThid),myByHi(myThid)
                0296          DO bi=myBxLo(myThid),myBxHi(myThid)
                0297           DO j=1-OLy,sNy+OLy
                0298            DO i=1-OLx,sNx+OLx
                0299             uIce(i,j,bi,bj) = uIce(i,j,bi,bj)*seaiceMaskU(i,j,bi,bj)
                0300             vIce(i,j,bi,bj) = vIce(i,j,bi,bj)*seaiceMaskV(i,j,bi,bj)
                0301            ENDDO
                0302           ENDDO
                0303          ENDDO
                0304         ENDDO
                0305 #endif /* SEAICE_CGRID */
                0306         CALL EXCH_UV_XY_RL( uIce, vIce, .TRUE., myThid )
                0307        ENDIF
                0308 
6060ec2938 Dimi*0309 C--   Read initial sea-ice thickness from file if available.
e7c33124ce Dimi*0310        IF ( HeffFile .NE. ' ' ) THEN
772590b63c Mart*0311         CALL READ_FLD_XY_RL( HeffFile, ' ', HEFF, 0, myThid )
                0312         _EXCH_XY_RL(HEFF,myThid)
e7c33124ce Dimi*0313         DO bj=myByLo(myThid),myByHi(myThid)
                0314          DO bi=myBxLo(myThid),myBxHi(myThid)
                0315           DO j=1-OLy,sNy+OLy
                0316            DO i=1-OLx,sNx+OLx
772590b63c Mart*0317             HEFF(i,j,bi,bj) = MAX(HEFF(i,j,bi,bj),ZERO)
e7c33124ce Dimi*0318            ENDDO
                0319           ENDDO
                0320          ENDDO
                0321         ENDDO
                0322        ENDIF
                0323 
                0324        DO bj=myByLo(myThid),myByHi(myThid)
                0325         DO bi=myBxLo(myThid),myBxHi(myThid)
                0326          DO j=1-OLy,sNy+OLy
                0327           DO i=1-OLx,sNx+OLx
772590b63c Mart*0328            IF(HEFF(i,j,bi,bj).GT.ZERO) AREA(i,j,bi,bj)=ONE
e7c33124ce Dimi*0329           ENDDO
                0330          ENDDO
                0331         ENDDO
                0332        ENDDO
f4cb8ca7db Jean*0333 
ff73878394 Dimi*0334 C--   Read initial sea-ice area from file if available.
09066b09cb Mart*0335        IF ( AreaFile .NE. ' ' ) THEN
772590b63c Mart*0336         CALL READ_FLD_XY_RL( AreaFile, ' ', AREA, 0, myThid )
                0337         _EXCH_XY_RL(AREA,myThid)
09066b09cb Mart*0338         DO bj=myByLo(myThid),myByHi(myThid)
                0339          DO bi=myBxLo(myThid),myBxHi(myThid)
                0340           DO j=1-OLy,sNy+OLy
                0341            DO i=1-OLx,sNx+OLx
772590b63c Mart*0342             AREA(i,j,bi,bj) = MAX(AREA(i,j,bi,bj),ZERO)
                0343             AREA(i,j,bi,bj) = MIN(AREA(i,j,bi,bj),ONE)
                0344             IF ( AREA(i,j,bi,bj) .LE. ZERO ) HEFF(i,j,bi,bj) = ZERO
                0345             IF ( HEFF(i,j,bi,bj) .LE. ZERO ) AREA(i,j,bi,bj) = ZERO
09066b09cb Mart*0346            ENDDO
                0347           ENDDO
                0348          ENDDO
                0349         ENDDO
                0350        ENDIF
                0351 
                0352        DO bj=myByLo(myThid),myByHi(myThid)
                0353         DO bi=myBxLo(myThid),myBxHi(myThid)
                0354          DO j=1-OLy,sNy+OLy
                0355           DO i=1-OLx,sNx+OLx
e5b8ebbabe Jean*0356            HSNOW(i,j,bi,bj) = 0.2 _d 0 * AREA(i,j,bi,bj)
09066b09cb Mart*0357           ENDDO
                0358          ENDDO
                0359         ENDDO
                0360        ENDDO
de31ea8481 Dimi*0361 
                0362 C--   Read initial snow thickness from file if available.
                0363        IF ( HsnowFile .NE. ' ' ) THEN
772590b63c Mart*0364         CALL READ_FLD_XY_RL( HsnowFile, ' ', HSNOW, 0, myThid )
                0365         _EXCH_XY_RL(HSNOW,myThid)
de31ea8481 Dimi*0366         DO bj=myByLo(myThid),myByHi(myThid)
                0367          DO bi=myBxLo(myThid),myBxHi(myThid)
                0368           DO j=1-OLy,sNy+OLy
                0369            DO i=1-OLx,sNx+OLx
772590b63c Mart*0370             HSNOW(i,j,bi,bj) = MAX(HSNOW(i,j,bi,bj),ZERO)
de31ea8481 Dimi*0371            ENDDO
                0372           ENDDO
                0373          ENDDO
                0374         ENDDO
                0375        ENDIF
fdfa8e151f Dimi*0376 
86b84a92fc Patr*0377 #ifdef SEAICE_ITD
                0378        DO bj=myByLo(myThid),myByHi(myThid)
                0379         DO bi=myBxLo(myThid),myBxHi(myThid)
                0380          DO j=1-OLy,sNy+OLy
                0381           DO i=1-OLx,sNx+OLx
5fe78992ba Mart*0382            AREAITD(i,j,1,bi,bj)   = AREA(i,j,bi,bj)
                0383            HEFFITD(i,j,1,bi,bj)   = HEFF(i,j,bi,bj)
                0384            HSNOWITD(i,j,1,bi,bj)  = HSNOW(i,j,bi,bj)
                0385            opnWtrFrac(i,j,bi,bj)  = 1. _d 0 - AREA(i,j,bi,bj)
                0386            fw2ObyRidge(i,j,bi,bj) = 0. _d 0
86b84a92fc Patr*0387           ENDDO
                0388          ENDDO
fd2abd71b8 Mart*0389          CALL SEAICE_ITD_REDIST(bi, bj, baseTime, nIter0, myThid)
86b84a92fc Patr*0390         ENDDO
                0391        ENDDO
                0392 #endif
                0393 
a98c4b8072 Ian *0394 #ifdef SEAICE_VARIABLE_SALINITY
3183247f03 Dimi*0395        DO bj=myByLo(myThid),myByHi(myThid)
                0396         DO bi=myBxLo(myThid),myBxHi(myThid)
f4cb8ca7db Jean*0397          DO j=1-OLy,sNy+OLy
                0398           DO i=1-OLx,sNx+OLx
0320e25227 Mart*0399            HSALT(i,j,bi,bj)=HEFF(i,j,bi,bj)*salt(i,j,kSrf,bi,bj)*
4eb4a54cba Jean*0400      &            SEAICE_rhoIce*SEAICE_saltFrac
                0401 cif     &            ICE2WATR*rhoConstFresh*SEAICE_saltFrac
a98c4b8072 Ian *0402 
3183247f03 Dimi*0403           ENDDO
                0404          ENDDO
                0405         ENDDO
                0406        ENDDO
                0407 
fdfa8e151f Dimi*0408 C--   Read initial sea ice salinity from file if available.
                0409        IF ( HsaltFile .NE. ' ' ) THEN
772590b63c Mart*0410         CALL READ_FLD_XY_RL( HsaltFile, ' ', HSALT, 0, myThid )
                0411         _EXCH_XY_RL(HSALT,myThid)
fdfa8e151f Dimi*0412        ENDIF
a98c4b8072 Ian *0413 #endif /* SEAICE_VARIABLE_SALINITY */
f4cb8ca7db Jean*0414 
e54fe3e1f9 Gael*0415 #ifdef ALLOW_SITRACER
f681b7f5d4 Dimi*0416 C--   Read initial sea ice age from file if available.
e54fe3e1f9 Gael*0417        DO iTr = 1, SItrMaxNum
                0418         IF ( SItrFile(iTr) .NE. ' ' ) THEN
                0419         CALL READ_FLD_XY_RL( siTrFile(iTr), ' ',
                0420      &   SItracer(1-OLx,1-OLy,1,1,iTr), 0, myThid )
                0421         _EXCH_XY_RL(SItracer(1-OLx,1-OLy,1,1,iTr),myThid)
ccaa3c61f4 Patr*0422         ENDIF
                0423        ENDDO
e54fe3e1f9 Gael*0424 #endif /* ALLOW_SITRACER */
f681b7f5d4 Dimi*0425 
809c36b928 Patr*0426       ENDIF
                0427 
4f6e464f38 Jean*0428 #ifdef ALLOW_OBCS
3352db77f4 Jean*0429 C--   In case we use scheme with a large stencil that extends into overlap:
                0430 C     no longer needed with the right masking in advection & diffusion S/R.
                0431 c     IF ( useOBCS ) THEN
                0432 c       DO bj=myByLo(myThid),myByHi(myThid)
                0433 c        DO bi=myBxLo(myThid),myBxHi(myThid)
48474b8e08 Jean*0434 c          CALL OBCS_COPY_TRACER( HEFF(1-OLx,1-OLy,bi,bj),
3352db77f4 Jean*0435 c    I                            1, bi, bj, myThid )
48474b8e08 Jean*0436 c          CALL OBCS_COPY_TRACER( AREA(1-OLx,1-OLy,bi,bj),
3352db77f4 Jean*0437 c    I                            1, bi, bj, myThid )
48474b8e08 Jean*0438 c          CALL OBCS_COPY_TRACER( HSNOW(1-OLx,1-OLy,bi,bj),
3352db77f4 Jean*0439 c    I                            1, bi, bj, myThid )
a98c4b8072 Ian *0440 #ifdef SEAICE_VARIABLE_SALINITY
48474b8e08 Jean*0441 c          CALL OBCS_COPY_TRACER( HSALT(1-OLx,1-OLy,bi,bj),
3352db77f4 Jean*0442 c    I                            1, bi, bj, myThid )
4f6e464f38 Jean*0443 #endif
3352db77f4 Jean*0444 c        ENDDO
                0445 c       ENDDO
                0446 c     ENDIF
4f6e464f38 Jean*0447 #endif /* ALLOW_OBCS */
                0448 
809c36b928 Patr*0449 C---  Complete initialization
                0450       DO bj=myByLo(myThid),myByHi(myThid)
                0451        DO bi=myBxLo(myThid),myBxHi(myThid)
45315406aa Mart*0452 #if ( defined SEAICE_CGRID || defined SEAICE_BGRID_DYNAMICS )
809c36b928 Patr*0453         DO j=1-OLy,sNy+OLy
                0454          DO i=1-OLx,sNx+OLx
8e32c48b8f Mart*0455           ZETA(i,j,bi,bj)        = HEFF(i,j,bi,bj)*(1.0 _d 11)
                0456           ETA(i,j,bi,bj)         = ZETA(i,j,bi,bj)/SEAICE_eccen**2
                0457           PRESS0(i,j,bi,bj)      = SEAICE_strength*HEFF(i,j,bi,bj)
ba6cfc5714 Mart*0458      &         *EXP(-SEAICE_cStar*(ONE-AREA(i,j,bi,bj)))
8e32c48b8f Mart*0459           SEAICE_zMax(i,j,bi,bj) = SEAICE_zetaMaxFac*PRESS0(i,j,bi,bj)
                0460           SEAICE_zMin(i,j,bi,bj) = SEAICE_zetaMin
                0461           PRESS0(i,j,bi,bj)      = PRESS0(i,j,bi,bj)*HEFFM(i,j,bi,bj)
809c36b928 Patr*0462          ENDDO
                0463         ENDDO
45315406aa Mart*0464 #endif
f3ce416a61 Jean*0465         IF ( useRealFreshWaterFlux .AND. .NOT.useThSIce ) THEN
                0466          DO j=1-OLy,sNy+OLy
                0467           DO i=1-OLx,sNx+OLx
772590b63c Mart*0468            sIceLoad(i,j,bi,bj) = HEFF(i,j,bi,bj)*SEAICE_rhoIce
13533b2b42 Mart*0469      &                         + HSNOW(i,j,bi,bj)*SEAICE_rhoSnow
f4cb8ca7db Jean*0470 
f3ce416a61 Jean*0471           ENDDO
                0472          ENDDO
                0473         ENDIF
809c36b928 Patr*0474        ENDDO
                0475       ENDDO
2f5e8addfd Mart*0476 #ifdef SEAICE_CGRID
                0477 C     compute tensile strength factor k: tensileStrength = k*PRESS
0dcda34642 Jean*0478 C     can be done in initialisation phase as long as it depends only
2f5e8addfd Mart*0479 C     on depth
                0480       IF ( SEAICE_tensilFac .NE. 0. _d 0 ) THEN
                0481        recip_tensilDepth = 0. _d 0
0dcda34642 Jean*0482        IF ( SEAICE_tensilDepth .GT. 0. _d 0 )
2f5e8addfd Mart*0483      &      recip_tensilDepth = 1. _d 0 / SEAICE_tensilDepth
                0484        DO bj=myByLo(myThid),myByHi(myThid)
                0485         DO bi=myBxLo(myThid),myBxHi(myThid)
                0486          DO j=1-OLy,sNy+OLy
                0487           DO i=1-OLx,sNx+OLx
                0488            tensileStrFac(i,j,bi,bj) = SEAICE_tensilFac*HEFFM(i,j,bi,bj)
                0489      &          *exp(-ABS(R_low(i,j,bi,bj))*recip_tensilDepth)
                0490           ENDDO
                0491          ENDDO
                0492         ENDDO
                0493        ENDDO
                0494       ENDIF
45315406aa Mart*0495 # if ( defined SEAICE_ALLOW_JFNK || defined SEAICE_ALLOW_KRYLOV )
                0496 C     Computing this metric cannot be done in S/R SEAICE_INIT_FIXED
                0497 C     where it belongs, because globalArea is only defined later after
                0498 C     S/R PACKAGES_INIT_FIXED, so we move this computation here.
                0499       CALL SEAICE_MAP_RS2VEC( nVec, rAw, rAs,
                0500      &     scalarProductMetric, .TRUE., myThid )
                0501       DO bj=myByLo(myThid),myByHi(myThid)
                0502        DO bi=myBxLo(myThid),myBxHi(myThid)
                0503         DO i=1,nVec
                0504          scalarProductMetric(i,1,bi,bj) =
                0505      &        scalarProductMetric(i,1,bi,bj)/globalArea
                0506         ENDDO
                0507        ENDDO
                0508       ENDDO
                0509 # endif /* SEAICE_ALLOW_JFNK or KRYLOV */
                0510 
2f5e8addfd Mart*0511 #endif /* SEAICE_CGRID */
809c36b928 Patr*0512 
c69ccc91fb antn*0513 #ifdef ALLOW_AUTODIFF
                0514 C-    end if useSEAICE block
                0515       ENDIF
                0516 #endif
                0517 
809c36b928 Patr*0518       RETURN
                0519       END