pkg/seaice/seaice_init_fixed.F

7e4f7741f6 Patr*0001 #include "SEAICE_OPTIONS.h"
                0002
5bb179ddc2 Mart*0003 CBOP
                0004 C !ROUTINE: SEAICE_INIT_FIXED
                0005
                0006 C !INTERFACE: ==========================================================
7e4f7741f6 Patr*0007       SUBROUTINE SEAICE_INIT_FIXED( myThid )
5bb179ddc2 Mart*0008 C !DESCRIPTION: \bv
afaeb4fe62 Jean*0009 C     *==========================================================*
                0010 C     | SUBROUTINE SEAICE_INIT_FIXED
                0011 C     | o Initialization of sea ice model.
                0012 C     *==========================================================*
                0013 C     *==========================================================*
7e4f7741f6 Patr*0014       IMPLICIT NONE
5bb179ddc2 Mart*0015 C \ev
afaeb4fe62 Jean*0016
5bb179ddc2 Mart*0017 C !USES: ===============================================================
7e4f7741f6 Patr*0018 #include "SIZE.h"
                0019 #include "EEPARAMS.h"
                0020 #include "PARAMS.h"
                0021 #include "GRID.h"
76d465095d Jean*0022 #include "FFIELDS.h"
03c669d1ab Jean*0023 #include "SEAICE_SIZE.h"
76d465095d Jean*0024 #include "SEAICE_PARAMS.h"
3f0f10fc37 Mart*0025 #include "SEAICE_GRID.h"
7e4f7741f6 Patr*0026 #include "SEAICE.h"
8bc8bee483 Gael*0027 #include "SEAICE_TRACER.h"
7e4f7741f6 Patr*0028
5bb179ddc2 Mart*0029 C !INPUT PARAMETERS: ===================================================
7e4f7741f6 Patr*0030 C     myThid - Thread no. that called this routine.
                0031       INTEGER myThid
                0032 CEndOfInterface
afaeb4fe62 Jean*0033
5bb179ddc2 Mart*0034 C !LOCAL VARIABLES: ====================================================
86b84a92fc Patr*0035 #ifdef SEAICE_ITD
                0036 C     msgBuf      :: Informational/error message buffer
                0037       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0038 C     k - loop counter for ITD categories
                0039       INTEGER k
514a01235c Jean*0040       _RL tmpVar
541c48397a Mart*0041       LOGICAL computeHlimit
86b84a92fc Patr*0042 #endif
ed6012c5a0 Jean*0043 C     i,j, bi,bj  :: Loop counters
ec0d7df165 Mart*0044       INTEGER i, j, bi, bj
0320e25227 Mart*0045       INTEGER kSrf
8bc8bee483 Gael*0046 #ifdef ALLOW_SITRACER
                0047       INTEGER iTracer
                0048 #endif
0320e25227 Mart*0049 C     local copy of surface layer thickness in meters
                0050       _RL dzSurf
45315406aa Mart*0051 #ifdef SEAICE_BGRID_DYNAMICS
5fe78992ba Mart*0052       _RL mask_uice
                0053 #endif
5bb179ddc2 Mart*0054 #ifdef SEAICE_ALLOW_SIDEDRAG
                0055       LOGICAL readCoastLineFields
                0056 #endif
                0057 CEOP
7e4f7741f6 Patr*0058
0320e25227 Mart*0059       IF ( usingPCoords ) THEN
                0060        kSrf        = Nr
84cb676f82 Mart*0061       ELSE
0320e25227 Mart*0062        kSrf        = 1
7e4f7741f6 Patr*0063       ENDIF
                0064
afaeb4fe62 Jean*0065 C     Initialize MNC variable information for SEAICE
b7411f1a84 Jean*0066       IF ( useMNC .AND. ( seaice_dump_mnc .OR. SEAICE_mon_mnc ) ) THEN
afaeb4fe62 Jean*0067         CALL SEAICE_MNC_INIT( myThid )
                0068       ENDIF
                0069
514a01235c Jean*0070 C     Only Master Thread updates parameter in commom block:
                0071       _BEGIN_MASTER(myThid)
                0072
6cbc659de0 Mart*0073 C     restart parameter
e501eee760 Mart*0074       SEAICEmomStartBDF = 0
                0075       IF ( SEAICEuseBDF2 ) SEAICEmomStartBDF = nIter0
6cbc659de0 Mart*0076
1c278edd09 Jean*0077 C--   Set mcPheePiston coeff (if still unset)
0320e25227 Mart*0078        dzSurf = drF(kSrf)
                0079        IF ( usingPCoords )
                0080      &      dzSurf = drF(kSrf) * recip_rhoConst * recip_gravity
1c278edd09 Jean*0081       IF ( SEAICE_mcPheePiston.EQ.UNSET_RL ) THEN
                0082         IF ( SEAICE_availHeatFrac.NE.UNSET_RL ) THEN
                0083           SEAICE_mcPheePiston = SEAICE_availHeatFrac
0320e25227 Mart*0084      &                        * dzSurf/SEAICE_deltaTtherm
1c278edd09 Jean*0085         ELSE
                0086           SEAICE_mcPheePiston = MCPHEE_TAPER_FAC
                0087      &                        * STANTON_NUMBER * USTAR_BASE
                0088           SEAICE_mcPheePiston = MIN( SEAICE_mcPheePiston,
0320e25227 Mart*0089      &                          dzSurf/SEAICE_deltaTtherm )
1c278edd09 Jean*0090         ENDIF
0320e25227 Mart*0091
1c278edd09 Jean*0092       ENDIF
                0093
8bc8bee483 Gael*0094 C--   SItracer specifications for basic tracers
                0095 #ifdef ALLOW_SITRACER
                0096       DO iTracer = 1, SItrNumInUse
1c278edd09 Jean*0097 C     "ice concentration" tracer that should remain .EQ.1.
                0098        IF (SItrName(iTracer).EQ.'one') THEN
8bc8bee483 Gael*0099          SItrFromOcean0(iTracer)    =ONE
                0100          SItrFromFlood0(iTracer)    =ONE
                0101          SItrExpand0(iTracer)       =ONE
                0102          SItrFromOceanFrac(iTracer) =ZERO
                0103          SItrFromFloodFrac(iTracer) =ZERO
1c278edd09 Jean*0104        ENDIF
                0105 C     age tracer: no age in ocean, or effect from ice cover changes
                0106        IF (SItrName(iTracer).EQ.'age') THEN
8bc8bee483 Gael*0107          SItrFromOcean0(iTracer)    =ZERO
                0108          SItrFromFlood0(iTracer)    =ZERO
                0109          SItrExpand0(iTracer)       =ZERO
                0110          SItrFromOceanFrac(iTracer) =ZERO
                0111          SItrFromFloodFrac(iTracer) =ZERO
1c278edd09 Jean*0112        ENDIf
                0113 C     salinity tracer:
                0114        IF (SItrName(iTracer).EQ.'salinity') THEN
8bc8bee483 Gael*0115          SItrMate(iTracer)          ='HEFF'
                0116          SItrExpand0(iTracer)       =ZERO
1c278edd09 Jean*0117          IF ( SEAICE_salinityTracer ) THEN
                0118            SEAICE_salt0    = ZERO
                0119            SEAICE_saltFrac = ZERO
                0120          ENDIF
                0121        ENDIF
                0122 C     simple, made up, ice surface roughness index prototype
                0123        IF (SItrName(iTracer).EQ.'ridge') THEN
8bc8bee483 Gael*0124          SItrMate(iTracer)          ='AREA'
                0125          SItrFromOcean0(iTracer)    =ZERO
                0126          SItrFromFlood0(iTracer)    =ZERO
                0127          SItrExpand0(iTracer)       =ZERO
                0128          SItrFromOceanFrac(iTracer) =ZERO
                0129          SItrFromFloodFrac(iTracer) =ZERO
1c278edd09 Jean*0130        ENDIF
f61838dfc1 Torg*0131 #ifdef SEAICE_GREASE
                0132 C     grease ice tracer:
                0133 c     (Smedrud and Martin, 2014, Ann. Glac.)
                0134        IF (SItrName(iTracer).EQ.'grease') THEN
                0135          SItrMate(iTracer)          ='HEFF'
                0136          SItrFromOcean0(iTracer)    =ZERO
                0137          SItrFromFlood0(iTracer)    =ZERO
                0138          SItrExpand0(iTracer)       =ZERO
                0139          SItrFromOceanFrac(iTracer) =ZERO
                0140          SItrFromFloodFrac(iTracer) =ZERO
                0141        ENDIF
                0142 #endif /* SEAICE_GREASE */
8bc8bee483 Gael*0143       ENDDO
f61838dfc1 Torg*0144 #endif /* ALLOW_SITRACER */
8bc8bee483 Gael*0145
86b84a92fc Patr*0146 #ifdef SEAICE_ITD
541c48397a Mart*0147 C     zeroth category needs to be zero
                0148       Hlimit(0)    = 0. _d 0
                0149 C     thickest category is unlimited
                0150       Hlimit(nITD) = 999.9 _d 0
                0151       computeHlimit=.FALSE.
                0152 C     check if Hlimit contains useful values
                0153       DO k = 1, nITD
                0154        IF ( Hlimit(k).EQ.UNSET_RL )    computeHlimit=.TRUE.
                0155        IF ( Hlimit(k).LE.Hlimit(k-1) ) computeHlimit=.TRUE.
                0156       ENDDO
                0157       IF ( computeHlimit ) THEN
22c740e9fe Mart*0158        WRITE(msgBuf,'(A,I2,A)') 'SEAICE_INIT_FIXED: Computing ',
541c48397a Mart*0159      &      nITD, ' thickness category limits with'
                0160        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0161      &                     SQUEEZE_RIGHT, myThid )
22c740e9fe Mart*0162        CALL WRITE_0D_RL( Hlimit_c1  ,INDEX_NONE,
                0163      &      'Hlimit_c1  =', ' /* ITD bin parameter */')
                0164        CALL WRITE_0D_RL( Hlimit_c2  ,INDEX_NONE,
                0165      &      'Hlimit_c2  =', ' /* ITD bin parameter */')
                0166        CALL WRITE_0D_RL( Hlimit_c3  ,INDEX_NONE,
                0167      &      'Hlimit_c3  =', ' /* ITD bin parameter */')
                0168        WRITE(msgBuf,'(A)') ' '
541c48397a Mart*0169        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0170      &                     SQUEEZE_RIGHT, myThid )
86b84a92fc Patr*0171 C     use Equ. 22 of Lipscomb et al. (2001, JGR) to generate ice
                0172 C     thickness category limits:
                0173 C     - dependends on given number of categories nITD
                0174 C     - choose between original parameters of Lipscomb et al. (2001):
becfda7837 Torg*0175 C       c1=3.0/N, c2=15*c1, c3=3.0
                0176 C       or emphasize thin end of ITD (in order to enhance ice growth):
                0177 C       c1=1.5/N, c2=42*c1, c3=3.3
                0178 C       -> HINT: set parameters c1, c2 and c3 in seaice_readparms.F
541c48397a Mart*0179        IF ( nITD.GT.1 ) THEN
                0180         tmpVar = nITD
                0181         tmpVar = oneRL / tmpVar
                0182         Hlimit_c1 = Hlimit_c1*tmpVar
                0183         Hlimit_c2 = Hlimit_c2*Hlimit_c1
                0184         DO k=1,nITD-1
                0185          Hlimit(k) = Hlimit(k-1)
                0186      &             + Hlimit_c1
                0187      &             + Hlimit_c2
                0188      &    *( oneRL + TANH( Hlimit_c3 *( FLOAT(k-1)*tmpVar - oneRL ) ) )
                0189         ENDDO
                0190        ENDIF
86b84a92fc Patr*0191       ENDIF
                0192 C     thickest category is unlimited
514a01235c Jean*0193       Hlimit(nITD) = 999.9 _d 0
86b84a92fc Patr*0194
514a01235c Jean*0195 #endif /* SEAICE_ITD */
86b84a92fc Patr*0196
514a01235c Jean*0197 C     Only Master Thread updates parameter in common block:
                0198       _END_MASTER(myThid)
84cb676f82 Mart*0199
ab972c4b6b Mart*0200 #ifdef SEAICE_ALLOW_JFNK
514a01235c Jean*0201 C     Only Master Thread updates parameter in common block:
                0202       _BEGIN_MASTER(myThid)
ab972c4b6b Mart*0203 C     initialise some diagnostic counters for the JFNK solver
                0204       totalNewtonIters   = 0
                0205       totalNewtonFails   = 0
                0206       totalKrylovIters   = 0
                0207       totalKrylovFails   = 0
                0208       totalJFNKtimeSteps = 0
514a01235c Jean*0209       _END_MASTER(myThid)
438b73e6d6 Mart*0210 C     this cannot be done here, because globalArea is only defined
                0211 C     after S/R PACKAGES_INIT_FIXED, so we move it to S/R SEAICE_INIT_VARIA
c363713d9d Jean*0212 CML      CALL SEAICE_MAP2VEC( nVec, rAw, rAs,
438b73e6d6 Mart*0213 CML     &     scalarProductMetric, .TRUE., myThid )
                0214 CML      DO bj=myByLo(myThid),myByHi(myThid)
                0215 CML       DO bi=myBxLo(myThid),myBxHi(myThid)
                0216 CML        DO i=1,nVec
c363713d9d Jean*0217 CML         scalarProductMetric(i,1,bi,bj) =
438b73e6d6 Mart*0218 CML     &        scalarProductMetric(i,1,bi,bj)/globalArea
                0219 CML        ENDDO
                0220 CML       ENDDO
                0221 CML      ENDDO
ab972c4b6b Mart*0222 #endif /* SEAICE_ALLOW_JFNK */
                0223
514a01235c Jean*0224 C--   all threads wait for master to finish initialisation of shared params
                0225       _BARRIER
                0226
84cb676f82 Mart*0227 #ifdef ALLOW_DIAGNOSTICS
                0228       IF ( useDiagnostics ) THEN
                0229         CALL SEAICE_DIAGNOSTICS_INIT( myThid )
                0230       ENDIF
                0231 #endif
                0232
1ed503f8a3 Gael*0233 C--   Summarise pkg/seaice configuration
                0234       CALL SEAICE_SUMMARY( myThid )
03c669d1ab Jean*0235
ec0d7df165 Mart*0236 C--   Initialise grid parameters that do no change during the integration
                0237
                0238 C--   Initialise grid info
                0239       DO bj=myByLo(myThid),myByHi(myThid)
                0240        DO bi=myBxLo(myThid),myBxHi(myThid)
148dd84005 jm-c 0241 C-    loops on tile indices bi,bj:
ec0d7df165 Mart*0242         DO j=1-OLy,sNy+OLy
                0243          DO i=1-OLx,sNx+OLx
0320e25227 Mart*0244           HEFFM  (i,j,bi,bj) = maskC(i,j,kSrf,bi,bj)
3f0f10fc37 Mart*0245 #ifdef SEAICE_CGRID
0320e25227 Mart*0246           SIMaskU(i,j,bi,bj) = maskW(i,j,kSrf,bi,bj)
                0247           SIMaskV(i,j,bi,bj) = maskS(i,j,kSrf,bi,bj)
3f0f10fc37 Mart*0248 #endif
ec0d7df165 Mart*0249          ENDDO
                0250         ENDDO
45315406aa Mart*0251 #ifdef SEAICE_BGRID_DYNAMICS
ec0d7df165 Mart*0252         DO j=1-OLy+1,sNy+OLy
                0253          DO i=1-OLx+1,sNx+OLx
                0254           UVM(i,j,bi,bj)=0. _d 0
                0255           mask_uice=HEFFM(i,j,  bi,bj)+HEFFM(i-1,j-1,bi,bj)
                0256      &             +HEFFM(i,j-1,bi,bj)+HEFFM(i-1,j,  bi,bj)
                0257           IF(mask_uice.GT.3.5 _d 0) UVM(i,j,bi,bj)=1. _d 0
                0258          ENDDO
                0259         ENDDO
45315406aa Mart*0260 #endif
ec0d7df165 Mart*0261
                0262 C     coefficients for metric terms
3f0f10fc37 Mart*0263 #if ( defined SEAICE_CGRID || defined SEAICE_BGRID_DYNAMICS )
ec0d7df165 Mart*0264         DO j=1-OLy,sNy+OLy
                0265          DO i=1-OLx,sNx+OLx
5fe78992ba Mart*0266           k1AtC(i,j,bi,bj) = 0.0 _d 0
                0267           k2AtC(i,j,bi,bj) = 0.0 _d 0
                0268           k1AtU(i,j,bi,bj) = 0.0 _d 0
                0269           k1AtV(i,j,bi,bj) = 0.0 _d 0
                0270           k2AtU(i,j,bi,bj) = 0.0 _d 0
                0271           k2AtV(i,j,bi,bj) = 0.0 _d 0
3f0f10fc37 Mart*0272 # ifdef SEAICE_CGRID
                0273           k1AtZ(i,j,bi,bj) = 0.0 _d 0
                0274           k2AtZ(i,j,bi,bj) = 0.0 _d 0
                0275 # endif
ec0d7df165 Mart*0276          ENDDO
                0277         ENDDO
3f0f10fc37 Mart*0278         IF ( SEAICEselectMetricTerms .GT. 0 ) THEN
                0279          IF ( usingSphericalPolarGrid ) THEN
ec0d7df165 Mart*0280 C     This is the only case where tan(phi) is not zero. In this case
                0281 C     C and U points, and Z and V points have the same phi, so that we
                0282 C     only need a copy here. Do not use tan(YC) and tan(YG), because
                0283 C     these
                0284 C     can be the geographical coordinates and not the correct grid
                0285 C     coordinates when the grid is rotated (phi/theta/psiEuler .NE. 0)
3f0f10fc37 Mart*0286           DO j=1-OLy,sNy+OLy
                0287            DO i=1-OLx,sNx+OLx
                0288             k2AtC(i,j,bi,bj) = - _tanPhiAtU(i,j,bi,bj)*recip_rSphere
                0289             k2AtU(i,j,bi,bj) = - _tanPhiAtU(i,j,bi,bj)*recip_rSphere
                0290             k2AtV(i,j,bi,bj) = - _tanPhiAtV(i,j,bi,bj)*recip_rSphere
                0291 # ifdef SEAICE_CGRID
                0292             k2AtZ(i,j,bi,bj) = - _tanPhiAtV(i,j,bi,bj)*recip_rSphere
                0293 # endif
                0294            ENDDO
ec0d7df165 Mart*0295           ENDDO
3f0f10fc37 Mart*0296          ELSEIF ( usingCurvilinearGrid ) THEN
ec0d7df165 Mart*0297 C     compute metric term coefficients from finite difference
                0298 C     approximation
3f0f10fc37 Mart*0299           DO j=1-OLy,sNy+OLy
                0300            DO i=1-OLx,sNx+OLx-1
                0301             k1AtC(i,j,bi,bj) = _recip_dyF(i,j,bi,bj)
                0302      &           * ( _dyG(i+1,j,bi,bj) - _dyG(i,j,bi,bj) )
                0303      &           * _recip_dxF(i,j,bi,bj)
                0304            ENDDO
ec0d7df165 Mart*0305           ENDDO
3f0f10fc37 Mart*0306           DO j=1-OLy,sNy+OLy-1
                0307            DO i=1-OLx,sNx+OLx
                0308             k2AtC(i,j,bi,bj) = _recip_dxF(i,j,bi,bj)
                0309      &           * ( _dxG(i,j+1,bi,bj) - _dxG(i,j,bi,bj) )
                0310      &           * _recip_dyF(i,j,bi,bj)
                0311            ENDDO
ec0d7df165 Mart*0312           ENDDO
3f0f10fc37 Mart*0313           DO j=1-OLy,sNy+OLy
                0314            DO i=1-OLx+1,sNx+OLx
                0315             k1AtU(i,j,bi,bj) = _recip_dyG(i,j,bi,bj)
                0316      &           * ( _dyF(i,j,bi,bj) - _dyF(i-1,j,bi,bj) )
                0317      &           * _recip_dxC(i,j,bi,bj)
                0318            ENDDO
ec0d7df165 Mart*0319           ENDDO
3f0f10fc37 Mart*0320           DO j=1-OLy,sNy+OLy-1
                0321            DO i=1-OLx,sNx+OLx
                0322             k2AtU(i,j,bi,bj) = _recip_dxC(i,j,bi,bj)
                0323      &           * ( _dxV(i,j+1,bi,bj) - _dxV(i,j,bi,bj) )
                0324      &           * _recip_dyG(i,j,bi,bj)
                0325            ENDDO
ec0d7df165 Mart*0326           ENDDO
3f0f10fc37 Mart*0327           DO j=1-OLy,sNy+OLy
                0328            DO i=1-OLx,sNx+OLx-1
                0329             k1AtV(i,j,bi,bj) = _recip_dyC(i,j,bi,bj)
                0330      &           * ( _dyU(i+1,j,bi,bj) - _dyU(i,j,bi,bj) )
                0331      &           * _recip_dxG(i,j,bi,bj)
                0332            ENDDO
ec0d7df165 Mart*0333           ENDDO
3f0f10fc37 Mart*0334           DO j=1-OLy+1,sNy+OLy
                0335            DO i=1-OLx,sNx+OLx
                0336             k2AtV(i,j,bi,bj) = _recip_dxG(i,j,bi,bj)
                0337      &           * ( _dxF(i,j,bi,bj) - _dxF(i,j-1,bi,bj) )
                0338      &           * _recip_dyC(i,j,bi,bj)
                0339            ENDDO
                0340           ENDDO
                0341 # ifdef SEAICE_CGRID
                0342           DO j=1-OLy,sNy+OLy
                0343            DO i=1-OLx+1,sNx+OLx
                0344             k1AtZ(i,j,bi,bj) = _recip_dyU(i,j,bi,bj)
                0345      &           * ( _dyC(i,j,bi,bj) - _dyC(i-1,j,bi,bj) )
                0346      &           * _recip_dxV(i,j,bi,bj)
                0347            ENDDO
                0348           ENDDO
                0349           DO j=1-OLy+1,sNy+OLy
                0350            DO i=1-OLx,sNx+OLx
                0351             k2AtZ(i,j,bi,bj) = _recip_dxV(i,j,bi,bj)
                0352      &           * ( _dxC(i,j,bi,bj) - _dxC(i,j-1,bi,bj) )
                0353      &           * _recip_dyU(i,j,bi,bj)
                0354            ENDDO
                0355           ENDDO
                0356 # endif /* SEAICE_CGRID */
                0357          ENDIF
                0358         ENDIF
                0359 #endif /* SEAICE_CGRID or SEAICE_BGRID_DYNAMICS */
                0360
                0361 #ifdef SEAICE_CGRID
                0362         IF ( SEAICEselectMetricTerms .LT. 2 ) THEN
                0363 C     Resetting metric factors at U,V-pointw to zero removes the extra
                0364 C     terms in the stress divergence
                0365          DO j=1-OLy,sNy+OLy
ec0d7df165 Mart*0366           DO i=1-OLx,sNx+OLx
3f0f10fc37 Mart*0367            k1AtU(i,j,bi,bj) = 0.0 _d 0
                0368            k2AtU(i,j,bi,bj) = 0.0 _d 0
                0369            k1AtV(i,j,bi,bj) = 0.0 _d 0
                0370            k2AtV(i,j,bi,bj) = 0.0 _d 0
ec0d7df165 Mart*0371           ENDDO
                0372          ENDDO
                0373         ENDIF
3f0f10fc37 Mart*0374 #endif /* SEAICE_CGRID */
                0375
                0376 #ifdef SEAICE_BGRID_DYNAMICS
45315406aa Mart*0377 C--   Choose a proxy level for geostrophic velocity,
                0378         DO j=1-OLy,sNy+OLy
                0379          DO i=1-OLx,sNx+OLx
                0380           KGEO(i,j,bi,bj)   = 0
                0381          ENDDO
                0382         ENDDO
                0383         DO j=1-OLy,sNy+OLy
                0384          DO i=1-OLx,sNx+OLx
                0385 #ifdef SEAICE_BICE_STRESS
                0386           KGEO(i,j,bi,bj) = 1
                0387 #else /* SEAICE_BICE_STRESS */
                0388           IF (klowc(i,j,bi,bj) .LT. 2) THEN
                0389            KGEO(i,j,bi,bj) = 1
                0390           ELSE
                0391            KGEO(i,j,bi,bj) = 2
                0392            DO WHILE ( abs(rC(KGEO(i,j,bi,bj))) .LT. 50.0 _d 0 .AND.
                0393      &          KGEO(i,j,bi,bj) .LT. (klowc(i,j,bi,bj)-1) )
                0394               KGEO(i,j,bi,bj) = KGEO(i,j,bi,bj) + 1
                0395            ENDDO
                0396           ENDIF
                0397 #endif /* SEAICE_BICE_STRESS */
                0398          ENDDO
                0399         ENDDO
                0400 #endif /* SEAICE_BGRID_DYNAMICS */
148dd84005 jm-c 0401 C-    end bi,bj loops
ec0d7df165 Mart*0402        ENDDO
                0403       ENDDO
                0404
5bb179ddc2 Mart*0405 #ifdef SEAICE_ALLOW_SIDEDRAG
                0406       DO bj=myByLo(myThid),myByHi(myThid)
                0407        DO bi=myBxLo(myThid),myBxHi(myThid)
                0408         DO j=1-OLy+1,sNy+OLy-1
                0409          DO i=1-OLx+1,sNx+OLx-1
                0410           coastRoughU(i,j,bi,bj) = ( 3. _d 0
                0411      &         - SIMaskU(i,j-1,bi,bj)
                0412      &         - SIMaskU(i,j,bi,bj)
                0413      &         - SIMaskU(i,j+1,bi,bj)) * SIMaskU(i,j,bi,bj)
                0414           coastRoughV(i,j,bi,bj) = ( 3. _d 0
                0415      &         - SIMaskV(i-1,j,bi,bj)
                0416      &         - SIMaskV(i,j,bi,bj)
                0417      &         - SIMaskV(i+1,j,bi,bj) ) * SIMaskV(i,j,bi,bj)
                0418          ENDDO
                0419         ENDDO
                0420        ENDDO
                0421       ENDDO
                0422
                0423       readCoastLineFields =
                0424      &     uCoastLineFile .NE. ' ' .OR. vCoastLineFile .NE. ' '
                0425       IF ( readCoastLineFields ) THEN
                0426        IF ( uCoastLineFile .NE. ' ' ) THEN
                0427         CALL READ_FLD_XY_RL( uCoastLineFile,' ',coastRoughU,0,myThid )
                0428         DO bj=myByLo(myThid),myByHi(myThid)
                0429          DO bi=myBxLo(myThid),myBxHi(myThid)
                0430           DO j=1-OLy,sNy+OLy-1
                0431            DO i=1-OLx,sNx+OLx
                0432             coastRoughU(i,j,bi,bj) = coastRoughU(i,j,bi,bj)
                0433      &           * SIMaskU(i,j,bi,bj) * recip_dxC(i,j,bi,bj)
                0434            ENDDO
                0435           ENDDO
                0436          ENDDO
                0437         ENDDO
                0438        ENDIF
                0439        IF ( vCoastLineFile .NE. ' ' ) THEN
                0440         CALL READ_FLD_XY_RL( vCoastLineFile,' ',coastRoughV,0,myThid )
                0441         DO bj=myByLo(myThid),myByHi(myThid)
                0442          DO bi=myBxLo(myThid),myBxHi(myThid)
                0443           DO j=1-OLy,sNy+OLy
                0444            DO i=1-OLx,sNx+OLx-1
                0445             coastRoughV(i,j,bi,bj) = coastRoughV(i,j,bi,bj)
                0446      &           * SIMaskV(i,j,bi,bj) * recip_dyC(i,j,bi,bj)
                0447            ENDDO
                0448           ENDDO
                0449          ENDDO
                0450         ENDDO
                0451        ENDIF
                0452       ENDIF
                0453
                0454       CALL EXCH_UV_XY_RL( coastRoughU, coastRoughV, .FALSE., myThid )
204068be6a R. S*0455 C     Write coastRoughU/V fields only we are actually going to use them
                0456 C     to reduce unnecessary output
                0457       IF ( SEAICEsideDrag .NE. 0. _d 0 ) THEN
                0458        CALL WRITE_FLD_XY_RL( 'coastRoughUnormalized',
                0459      &      ' ', coastRoughU, -1, myThid )
                0460        CALL WRITE_FLD_XY_RL( 'coastRoughVnormalized',
                0461      &      ' ', coastRoughV, -1, myThid )
                0462       ENDIF
5bb179ddc2 Mart*0463 #endif /* SEAICE_ALLOW_SIDEDRAG */
                0464
148dd84005 jm-c 0465 #ifdef ALLOW_SHELFICE
                0466       IF ( useShelfIce ) THEN
                0467        DO bj=myByLo(myThid),myByHi(myThid)
                0468         DO bi=myBxLo(myThid),myBxHi(myThid)
                0469 C--   Prevent seaice to form where Ice-Shelf is present
                0470           CALL SHELFICE_MASK_SEAICE(
                0471      U                  HEFFM,
                0472      I                  bi, bj, -1, myThid )
                0473           DO j=2-OLy,sNy+OLy
                0474            DO i=2-OLx,sNx+OLx
                0475             IF ( HEFFM(i-1,j,bi,bj).EQ.zeroRL .OR.
                0476      &           HEFFM( i, j,bi,bj).EQ.zeroRL )
                0477      &           SIMaskU(i,j,bi,bj) = zeroRL
                0478             IF ( HEFFM(i,j-1,bi,bj).EQ.zeroRL .OR.
                0479      &           HEFFM(i, j, bi,bj).EQ.zeroRL )
                0480      &           SIMaskV(i,j,bi,bj) = zeroRL
                0481            ENDDO
                0482           ENDDO
                0483         ENDDO
                0484        ENDDO
                0485        CALL EXCH_UV_XY_RL( SIMaskU, SIMaskV, .FALSE., myThid )
                0486       ENDIF
                0487 #endif /* ALLOW_SHELFICE */
                0488
7e4f7741f6 Patr*0489       RETURN
                0490       END