Back to home page

MITgcm
 
 

    

File indexing completed on 2023-08-04 05:10:46 UTC

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

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!