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