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