File indexing completed on 2024-10-18 05:11:26 UTC
view on githubraw file Latest commit 5bb179dd on 2024-10-17 18:00:27 UTC
53092bcb42 Mart*0001 #include "SEAICE_OPTIONS.h"
fec75090d3 Jean*0002 #ifdef ALLOW_OBCS
0003 # include "OBCS_OPTIONS.h"
0004 #else
0005 # define OBCS_UVICE_OLD
0006 #endif
772b2ed80e Gael*0007 #ifdef ALLOW_AUTODIFF
0008 # include "AUTODIFF_OPTIONS.h"
0009 #endif
53092bcb42 Mart*0010
efcb5efb58 Jean*0011
0012
0013
0014
5fb12930f9 Mart*0015
bfe6f7447e Mart*0016
0017
5fb12930f9 Mart*0018
0019
efcb5efb58 Jean*0020
86fa82a64d Jean*0021
e7ae128053 Jean*0022
86fa82a64d Jean*0023
bb3037b865 Gael*0024 SUBROUTINE SEAICE_LSR( myTime, myIter, myThid )
86fa82a64d Jean*0025
0026
0027
0028
0029
0030
0031
0032
0033
efcb5efb58 Jean*0034
86fa82a64d Jean*0035
0036
0037
efcb5efb58 Jean*0038
86fa82a64d Jean*0039
0040
0041
589eca759f Mart*0042 IMPLICIT NONE
0043
0044
0045 #include "SIZE.h"
0046 #include "EEPARAMS.h"
0047 #include "PARAMS.h"
0048 #include "DYNVARS.h"
0049 #include "GRID.h"
9b4636ff01 Patr*0050 #include "SEAICE_SIZE.h"
589eca759f Mart*0051 #include "SEAICE_PARAMS.h"
03c669d1ab Jean*0052 #include "SEAICE.h"
589eca759f Mart*0053
0054 #ifdef ALLOW_AUTODIFF_TAMC
0055 # include "tamc.h"
0056 #endif
0057
86fa82a64d Jean*0058
589eca759f Mart*0059
86fa82a64d Jean*0060
0061
0062
589eca759f Mart*0063 _RL myTime
0064 INTEGER myIter
0065 INTEGER myThid
86fa82a64d Jean*0066
589eca759f Mart*0067
0068 #ifdef SEAICE_CGRID
7b7a25aa58 Jean*0069
0070 LOGICAL DIFFERENT_MULTIPLE
0071 EXTERNAL DIFFERENT_MULTIPLE
589eca759f Mart*0072
86fa82a64d Jean*0073
589eca759f Mart*0074
86fa82a64d Jean*0075
589eca759f Mart*0076
4a08d54d3a Mart*0077 INTEGER ipass
5fb12930f9 Mart*0078 INTEGER i, j, m, bi, bj
0079 INTEGER iMin, iMax, jMin, jMax
79df32c3f1 Mart*0080 INTEGER ICOUNT1, ICOUNT2, linearIterLoc, nonLinIterLoc
486cba1d02 Mart*0081 INTEGER kSrf
7b7a25aa58 Jean*0082 LOGICAL doIterate4u, doIterate4v
143d9ce879 Dami*0083 LOGICAL doNonLinLoop
86fa82a64d Jean*0084 CHARACTER*(MAX_LEN_MBUF) msgBuf
589eca759f Mart*0085
143d9ce879 Dami*0086 #ifdef ALLOW_DIAGNOSTICS
0087 LOGICAL DIAGNOSTICS_IS_ON
0088 EXTERNAL DIAGNOSTICS_IS_ON
0089 _RL RTmp (1:sNx,1:sNy)
0090 #endif /* ALLOW_DIAGNOSTICS */
0091
e45202e340 Mart*0092 _RL WFAU, WFAV, WFAU2, WFAV2
5fb12930f9 Mart*0093 _RL S1, S2, S1A, S2A
e232688b80 Mart*0094 _RL recip_deltaT
589eca759f Mart*0095
143d9ce879 Dami*0096 #ifdef SEAICE_ALLOW_LSR_FLEX
0097 _RL residIni, residIniNonLin
0098 _RL LSRflexFac, LSR_ERRORfu, LSR_ERRORfv
0099 _RL residkm1, residkm2
0100 #endif /* SEAICE_ALLOW_LSR_FLEX */
0101
589eca759f Mart*0102
03c669d1ab Jean*0103 _RL AU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0104 _RL BU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0105 _RL CU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0106 _RL AV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0107 _RL BV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0108 _RL CV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
2d6335c1f7 Jean*0109
03c669d1ab Jean*0110 _RL rhsU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0111 _RL rhsV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
2d6335c1f7 Jean*0112 #ifdef SEAICE_GLOBAL_3DIAG_SOLVER
0113
0114 _RL rhsUx(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0115 _RL rhsVy(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0116 INTEGER iSubIter
b2e217c119 Jean*0117 INTEGER errCode
2d6335c1f7 Jean*0118 #endif
7b7a25aa58 Jean*0119
03c669d1ab Jean*0120 _RL uRt1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0121 _RL uRt2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
7b7a25aa58 Jean*0122
03c669d1ab Jean*0123 _RL vRt1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0124 _RL vRt2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
589eca759f Mart*0125
03c669d1ab Jean*0126 _RL etaPlusZeta (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0127 _RL zetaMinusEta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
df1dac8b7b Mart*0128
0129 _RL dragSym (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
bfe6f7447e Mart*0130
0131 _RL uIceC(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0132 _RL vIceC(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
589eca759f Mart*0133
03c669d1ab Jean*0134 _RL uTmp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0135 _RL vTmp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
4a08d54d3a Mart*0136
0137 _RL fxTmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0138 _RL fyTmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
70e078b38a Mart*0139
733ca03edd Mart*0140 _RL areaW(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0141 _RL areaS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
8ce59fd29e Mart*0142 #ifdef SEAICE_ALLOW_MOM_ADVECTION
0143
0144 _RL gUmom(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
0145 _RL gVmom(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
0146 #endif /* SEAICE_ALLOW_MOM_ADVECTION */
d24acf3cfc Jean*0147 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0148
0149
0150
0151 INTEGER nlkey, lnkey, tkey
d24acf3cfc Jean*0152 #endif
35fda33b05 Jean*0153 _RL COSWAT
0154 _RS SINWAT
589eca759f Mart*0155 _RL UERR
63a26f3b3e Mart*0156 #ifdef SEAICE_ALLOW_CHECK_LSR_CONVERGENCE
0157 _RL resnorm, EKnorm, counter
0158 #endif /* SEAICE_ALLOW_CHECK_LSR_CONVERGENCE */
efcb5efb58 Jean*0159 LOGICAL printResidual
0160 _RL residUini, residVini, residUend, residVend
0161 #ifdef SEAICE_ALLOW_FREEDRIFT
03c669d1ab Jean*0162 _RL uRes (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0163 _RL vRes (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
efcb5efb58 Jean*0164 _RL errIni, errFD, errSum
0165 _RL residU_fd, residV_fd
0166 _RL residUmix, residVmix
0167 #endif /* SEAICE_ALLOW_FREEDRIFT */
0c560a7486 Mart*0168 _RL areaMinLoc
efcb5efb58 Jean*0169
0170
0171
0172 printResidual = debugLevel.GE.debLevA
0173 & .AND. DIFFERENT_MULTIPLE( SEAICE_monFreq, myTime, deltaTClock )
589eca759f Mart*0174
a68248c150 Mart*0175
0176 jMin = 1 - SEAICE_OLy
0177 jMax = sNy + SEAICE_OLy
0178 iMin = 1 - SEAICE_OLx
0179 iMax = sNx + SEAICE_OLx
cbf501ab81 Jean*0180
79df32c3f1 Mart*0181 linearIterLoc = SEAICElinearIterMax
0182 nonLinIterLoc = SEAICEnonLinIterMax
0183 IF ( SEAICEusePicardAsPrecon ) THEN
0184 linearIterLoc = SEAICEpreconlinIter
0185 nonLinIterLoc = SEAICEpreconNL_Iter
0186 ENDIF
143d9ce879 Dami*0187
e232688b80 Mart*0188 recip_deltaT = 1. _d 0 / SEAICE_deltaTdyn
5fb12930f9 Mart*0189
143d9ce879 Dami*0190
0191 IF ( usingPCoords ) THEN
0192 kSrf = Nr
0193 ELSE
0194 kSrf = 1
0195 ENDIF
0196
0197 SINWAT=SIN(SEAICE_waterTurnAngle*deg2rad)
0198 COSWAT=COS(SEAICE_waterTurnAngle*deg2rad)
0199
8377b8ee87 Mart*0200 #ifdef ALLOW_AUTODIFF
d976cd8bb3 Patr*0201
0202 DO bj=myByLo(myThid),myByHi(myThid)
0203 DO bi=myBxLo(myThid),myBxHi(myThid)
0204 DO j=1-OLy,sNy+OLy
0205 DO i=1-OLx,sNx+OLx
a16f976091 Mart*0206 deltaC (i,j,bi,bj) = 0. _d 0
0207 press (i,j,bi,bj) = 0. _d 0
7c0da9bbfa Mart*0208 zeta (i,j,bi,bj) = 0. _d 0
0fd6b2de1a Mart*0209 zetaZ (i,j,bi,bj) = 0. _d 0
7c0da9bbfa Mart*0210 eta (i,j,bi,bj) = 0. _d 0
0211 etaZ (i,j,bi,bj) = 0. _d 0
0212 uIceC (i,j,bi,bj) = 0. _d 0
0213 vIceC (i,j,bi,bj) = 0. _d 0
0214 uIceNm1(i,j,bi,bj) = 0. _d 0
0215 vIceNm1(i,j,bi,bj) = 0. _d 0
d976cd8bb3 Patr*0216 ENDDO
0217 ENDDO
0218 ENDDO
0219 ENDDO
0220 #endif
0221
70e078b38a Mart*0222
733ca03edd Mart*0223 DO bj=myByLo(myThid),myByHi(myThid)
0224 DO bi=myBxLo(myThid),myBxHi(myThid)
0225 DO j=1-OLy,sNy+OLy
0226 DO i=1-OLx,sNx+OLx
4a08d54d3a Mart*0227 areaW(i,j,bi,bj) = 1. _d 0
0228 areaS(i,j,bi,bj) = 1. _d 0
0229 ENDDO
733ca03edd Mart*0230 ENDDO
0c560a7486 Mart*0231 areaMinLoc = 0. _d 0
733ca03edd Mart*0232 IF ( SEAICEscaleSurfStress ) THEN
0c560a7486 Mart*0233 areaMinLoc = 1. _d -10
0234
4a08d54d3a Mart*0235 DO j=jMin,jMax
0236 DO i=iMin,iMax
0237 areaW(i,j,bi,bj) =
0238 & 0.5 _d 0*(AREA(i,j,bi,bj)+AREA(i-1,j,bi,bj))
0239 areaS(i,j,bi,bj) =
0240 & 0.5 _d 0*(AREA(i,j,bi,bj)+AREA(i,j-1,bi,bj))
733ca03edd Mart*0241 ENDDO
0242 ENDDO
0243 ENDIF
70e078b38a Mart*0244 ENDDO
0245 ENDDO
0246
4a08d54d3a Mart*0247 #ifdef ALLOW_AUTODIFF_TAMC
0248
0249
0250 #endif /* ALLOW_AUTODIFF_TAMC */
0251
0252
0253
0254
0255 DO bj=myByLo(myThid),myByHi(myThid)
0256 DO bi=myBxLo(myThid),myBxHi(myThid)
0257 DO j=1-OLy,sNy+OLy
0258 DO i=1-OLx,sNx+OLx
0259 uIceNm1(i,j,bi,bj) = uIce(i,j,bi,bj)
0260 vIceNm1(i,j,bi,bj) = vIce(i,j,bi,bj)
0261 fxTmp (i,j,bi,bj) = FORCEX0(i,j,bi,bj)
0262 & +seaiceMassU(i,j,bi,bj)*recip_deltaT
0263 & *uIceNm1(i,j,bi,bj)
0264 fyTmp (i,j,bi,bj) = FORCEY0(i,j,bi,bj)
0265 & +seaiceMassV(i,j,bi,bj)*recip_deltaT
0266 & *vIceNm1(i,j,bi,bj)
0267 ENDDO
0268 ENDDO
0269 ENDDO
0270 ENDDO
0d75a51072 Mart*0271
143d9ce879 Dami*0272 #ifdef SEAICE_ALLOW_LSR_FLEX
0273 IF ( SEAICEuseLSRflex ) THEN
0274 residkm1 = 1. _d 0
0275 residkm2 = 1. _d 0
0276 residIniNonLin = 1. _d 0
0277 ENDIF
0278 #endif /* SEAICE_ALLOW_LSR_FLEX */
0279
0280 doNonLinLoop = .TRUE.
0281
0282
0283
0284
0285
0286
0287
79df32c3f1 Mart*0288 DO ipass=1,nonLinIterLoc
9b4636ff01 Patr*0289
0290 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0291 nlkey = (ikey_dynamics-1)*MPSEUDOTIMESTEPS + ipass
9b4636ff01 Patr*0292 # ifdef SEAICE_ALLOW_FREEDRIFT
edb6656069 Mart*0293
0294
9b4636ff01 Patr*0295 # endif
edb6656069 Mart*0296
9b4636ff01 Patr*0297 #endif /* ALLOW_AUTODIFF_TAMC */
0298
143d9ce879 Dami*0299 IF ( doNonLinLoop ) THEN
bb3037b865 Gael*0300 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0301
0302
0303
0304
bb3037b865 Gael*0305 #endif /* ALLOW_AUTODIFF_TAMC */
0306
589eca759f Mart*0307
e45202e340 Mart*0308 WFAU=SEAICE_LSRrelaxU
0309 WFAV=SEAICE_LSRrelaxV
589eca759f Mart*0310 WFAU2=ZERO
0311 WFAV2=ZERO
7b7a25aa58 Jean*0312 S1 = 0.
0313 S2 = 0.
e45202e340 Mart*0314 S1A=0.80 _d 0
0315 S2A=0.80 _d 0
589eca759f Mart*0316
4a08d54d3a Mart*0317 DO bj=myByLo(myThid),myByHi(myThid)
0318 DO bi=myBxLo(myThid),myBxHi(myThid)
0319 IF ( ipass .EQ. 1 ) THEN
0320
a10601a61f Mart*0321 DO j=1-OLy,sNy+OLy
0322 DO i=1-OLx,sNx+OLx
4a08d54d3a Mart*0323 uIceC(i,j,bi,bj) = uIce(i,j,bi,bj)
0324 vIceC(i,j,bi,bj) = vIce(i,j,bi,bj)
a10601a61f Mart*0325 ENDDO
0326 ENDDO
4a08d54d3a Mart*0327 ELSEIF ( ipass .EQ. 2 .AND. SEAICEnonLinIterMax .LE. 2 ) THEN
0328
a10601a61f Mart*0329 DO j=1-OLy,sNy+OLy
0330 DO i=1-OLx,sNx+OLx
4a08d54d3a Mart*0331 uIce (i,j,bi,bj)=HALF*( uIce(i,j,bi,bj)+uIceNm1(i,j,bi,bj) )
0332 vIce (i,j,bi,bj)=HALF*( vIce(i,j,bi,bj)+vIceNm1(i,j,bi,bj) )
0333 uIceC(i,j,bi,bj)=uIce(i,j,bi,bj)
0334 vIceC(i,j,bi,bj)=vIce(i,j,bi,bj)
a10601a61f Mart*0335 ENDDO
0336 ENDDO
4a08d54d3a Mart*0337 ELSE
0338
0339
0340
de36f9bd7d Mart*0341 DO j=1-OLy,sNy+OLy
0342 DO i=1-OLx,sNx+OLx
4a08d54d3a Mart*0343
0344
0345
0346
0347 uIceC(i,j,bi,bj) = HALF*( uIce(i,j,bi,bj)+uIceC(i,j,bi,bj) )
0348 vIceC(i,j,bi,bj) = HALF*( vIce(i,j,bi,bj)+vIceC(i,j,bi,bj) )
de36f9bd7d Mart*0349 ENDDO
0350 ENDDO
4a08d54d3a Mart*0351 ENDIF
0352
de36f9bd7d Mart*0353 ENDDO
4a08d54d3a Mart*0354 ENDDO
a10601a61f Mart*0355
589eca759f Mart*0356 #ifdef ALLOW_AUTODIFF_TAMC
3e8bf7e743 Mart*0357
143d9ce879 Dami*0358
0359
3e8bf7e743 Mart*0360
edb6656069 Mart*0361
0362
0363
0364
589eca759f Mart*0365 #endif /* ALLOW_AUTODIFF_TAMC */
0366
0367 CALL SEAICE_CALC_STRAINRATES(
0368 I uIceC, vIceC,
0369 O e11, e22, e12,
3e8bf7e743 Mart*0370 I ipass, myTime, myIter, myThid )
589eca759f Mart*0371
0372 CALL SEAICE_CALC_VISCOSITIES(
8e32c48b8f Mart*0373 I e11, e22, e12, SEAICE_zMin, SEAICE_zMax, HEFFM, press0,
0374 I tensileStrFac,
0fd6b2de1a Mart*0375 O eta, etaZ, zeta, zetaZ, press, deltaC,
3e8bf7e743 Mart*0376 I ipass, myTime, myIter, myThid )
589eca759f Mart*0377
3f31c7d5de Mart*0378 CALL SEAICE_OCEANDRAG_COEFFS(
ec0d7df165 Mart*0379 I uIceC, vIceC, HEFFM,
3f31c7d5de Mart*0380 O DWATN,
3e8bf7e743 Mart*0381 I ipass, myTime, myIter, myThid )
df1dac8b7b Mart*0382
0383 #ifdef SEAICE_ALLOW_BOTTOMDRAG
d5254b4e3d Mart*0384 CALL SEAICE_BOTTOMDRAG_COEFFS(
ec0d7df165 Mart*0385 I uIceC, vIceC, HEFFM,
df1dac8b7b Mart*0386 #ifdef SEAICE_ITD
0387 I HEFFITD, AREAITD, AREA,
0388 #else
0389 I HEFF, AREA,
cbf501ab81 Jean*0390 #endif
df1dac8b7b Mart*0391 O CbotC,
0392 I ipass, myTime, myIter, myThid )
0393 #endif /* SEAICE_ALLOW_BOTTOMDRAG */
5bb179ddc2 Mart*0394
0395 #ifdef SEAICE_ALLOW_SIDEDRAG
0396 IF ( SEAICEsideDrag .NE. 0. _d 0 ) CALL SEAICE_SIDEDRAG_STRESS(
0397 I uIceC, vIceC, coastRoughU, coastRoughV, AREA,
0398 O sideDragU, sideDragV,
0399 I ipass, myTime, myIter, myThid )
0400 #ifdef ALLOW_AUTODIFF_TAMC
0401
0402
0403 #endif /* ALLOW_AUTODIFF_TAMC */
0404 #endif /* SEAICE_ALLOW_SIDEDRAG */
d681af2948 Mart*0405
0406
0407
0408 DO bj=myByLo(myThid),myByHi(myThid)
0409 DO bi=myBxLo(myThid),myBxHi(myThid)
4a08d54d3a Mart*0410 DO j=jMin-1,jMax
0411 DO i=iMin-1,iMax
0412 etaPlusZeta (i,j,bi,bj) = ETA (i,j,bi,bj)+ZETA(i,j,bi,bj)
0413 zetaMinusEta(i,j,bi,bj) = ZETA(i,j,bi,bj)-ETA (i,j,bi,bj)
d681af2948 Mart*0414 ENDDO
0415 ENDDO
4a08d54d3a Mart*0416 DO j=1-OLy,sNy+OLy
0417 DO i=1-OLx,sNx+OLx
0418 dragSym(i,j,bi,bj) = DWATN(i,j,bi,bj)*COSWAT
df1dac8b7b Mart*0419 #ifdef SEAICE_ALLOW_BOTTOMDRAG
4a08d54d3a Mart*0420 & +CbotC(i,j,bi,bj)
df1dac8b7b Mart*0421 #endif /* SEAICE_ALLOW_BOTTOMDRAG */
0422 ENDDO
0423 ENDDO
d681af2948 Mart*0424 ENDDO
0425 ENDDO
0426
0427
0428
0429
589eca759f Mart*0430 DO bj=myByLo(myThid),myByHi(myThid)
0431 DO bi=myBxLo(myThid),myBxHi(myThid)
4a08d54d3a Mart*0432 DO j=jMin,jMax
0433 DO i=iMin,iMax
efcb5efb58 Jean*0434
589eca759f Mart*0435
4a08d54d3a Mart*0436 FORCEX(i,j,bi,bj) = fxTmp(i,j,bi,bj)+
0437 & ( 0.5 _d 0 * ( DWATN(i,j,bi,bj)+DWATN(i-1,j,bi,bj) ) *
0438 & COSWAT * uVel(i,j,kSrf,bi,bj)
0439 & - SIGN(SINWAT, _fCori(i,j,bi,bj))* 0.5 _d 0 *
0440 & ( DWATN(i ,j,bi,bj) * 0.5 _d 0 *
0441 & (vVel(i ,j ,kSrf,bi,bj)-vIceC(i ,j ,bi,bj)
0442 & +vVel(i ,j+1,kSrf,bi,bj)-vIceC(i ,j+1,bi,bj))
0443 & + DWATN(i-1,j,bi,bj) * 0.5 _d 0 *
0444 & (vVel(i-1,j ,kSrf,bi,bj)-vIceC(i-1,j ,bi,bj)
0445 & +vVel(i-1,j+1,kSrf,bi,bj)-vIceC(i-1,j+1,bi,bj))
0446 & ) ) * areaW(i,j,bi,bj)
0447 FORCEY(i,j,bi,bj) = fyTmp(i,j,bi,bj)+
0448 & ( 0.5 _d 0 * ( DWATN(i,j,bi,bj)+DWATN(i,j-1,bi,bj) ) *
0449 & COSWAT * vVel(i,j,kSrf,bi,bj)
0450 & + SIGN(SINWAT, _fCori(i,j,bi,bj)) * 0.5 _d 0 *
0451 & ( DWATN(i,j ,bi,bj) * 0.5 _d 0 *
0452 & (uVel(i ,j ,kSrf,bi,bj)-uIceC(i ,j ,bi,bj)
0453 & +uVel(i+1,j ,kSrf,bi,bj)-uIceC(i+1,j ,bi,bj))
0454 & + DWATN(i,j-1,bi,bj) * 0.5 _d 0 *
0455 & (uVel(i ,j-1,kSrf,bi,bj)-uIceC(i ,j-1,bi,bj)
0456 & +uVel(i+1,j-1,kSrf,bi,bj)-uIceC(i+1,j-1,bi,bj))
0457 & ) ) * areaS(i,j,bi,bj)
589eca759f Mart*0458 ENDDO
0459 ENDDO
5fb12930f9 Mart*0460 DO j=jMin,jMax
0461 DO i=iMin,iMax
efcb5efb58 Jean*0462
4a08d54d3a Mart*0463 FORCEX(i,j,bi,bj) = FORCEX(i,j,bi,bj) + HALF*
0464 & ( seaiceMassC(i ,j,bi,bj) * _fCori(i ,j,bi,bj)
efcb5efb58 Jean*0465 & *0.5 _d 0*( vIceC( i ,j,bi,bj)+vIceC( i ,j+1,bi,bj) )
4a08d54d3a Mart*0466 & + seaiceMassC(i-1,j,bi,bj) * _fCori(i-1,j,bi,bj)
efcb5efb58 Jean*0467 & *0.5 _d 0*( vIceC(i-1,j,bi,bj)+vIceC(i-1,j+1,bi,bj) ) )
4a08d54d3a Mart*0468 FORCEY(i,j,bi,bj) = FORCEY(i,j,bi,bj) - HALF*
0469 & ( seaiceMassC(i,j ,bi,bj) * _fCori(i,j ,bi,bj)
efcb5efb58 Jean*0470 & *0.5 _d 0*( uIceC(i ,j ,bi,bj)+uIceC(i+1, j,bi,bj) )
4a08d54d3a Mart*0471 & + seaiceMassC(i,j-1,bi,bj) * _fCori(i,j-1,bi,bj)
efcb5efb58 Jean*0472 & *0.5 _d 0*( uIceC(i ,j-1,bi,bj)+uIceC(i+1,j-1,bi,bj) ) )
0473 ENDDO
0474 ENDDO
4a08d54d3a Mart*0475
5fb12930f9 Mart*0476 DO j=jMin,jMax
0477 DO i=iMin,iMax
4a08d54d3a Mart*0478 FORCEX(i,j,bi,bj) = FORCEX(i,j,bi,bj)*seaiceMaskU(i,j,bi,bj)
0479 FORCEY(i,j,bi,bj) = FORCEY(i,j,bi,bj)*seaiceMaskV(i,j,bi,bj)
589eca759f Mart*0480 ENDDO
0481 ENDDO
7b7a25aa58 Jean*0482
b7c2bb63b7 Mart*0483
0484
0485
bfe6f7447e Mart*0486 DO j=jMin,jMax
0487 DO i=iMin,iMax
4a08d54d3a Mart*0488 rhsU(i,j,bi,bj) = FORCEX(i,j,bi,bj)
3ff8c9e192 Jean*0489 ENDDO
0490 ENDDO
bfe6f7447e Mart*0491 CALL SEAICE_LSR_RHSU(
0fd6b2de1a Mart*0492 I zetaMinusEta, etaPlusZeta, etaZ, zetaZ, press,
cbf501ab81 Jean*0493 I uIceC, vIceC,
bfe6f7447e Mart*0494 U rhsU,
0495 I iMin, iMax, jMin, jMax, bi, bj, myThid )
b7c2bb63b7 Mart*0496
0497
0498
bfe6f7447e Mart*0499 DO j=jMin,jMax
0500 DO i=iMin,iMax
4a08d54d3a Mart*0501 rhsV(i,j,bi,bj) = FORCEY(i,j,bi,bj)
3ff8c9e192 Jean*0502 ENDDO
0503 ENDDO
bfe6f7447e Mart*0504 CALL SEAICE_LSR_RHSV(
0fd6b2de1a Mart*0505 I zetaMinusEta, etaPlusZeta, etaZ, zetaZ, press,
cbf501ab81 Jean*0506 I uIceC, vIceC,
bfe6f7447e Mart*0507 U rhsV,
0508 I iMin, iMax, jMin, jMax, bi, bj, myThid )
8ce59fd29e Mart*0509
0510 #ifdef SEAICE_ALLOW_MOM_ADVECTION
cbf501ab81 Jean*0511 IF ( SEAICEmomAdvection ) THEN
4a08d54d3a Mart*0512 DO j=1-OLy,sNy+OLy
0513 DO i=1-OLx,sNx+OLx
0514 gUmom(i,j) = 0. _d 0
0515 gVmom(i,j) = 0. _d 0
8ce59fd29e Mart*0516 ENDDO
0517 ENDDO
0518 CALL SEAICE_MOM_ADVECTION(
0519 I bi,bj,iMin,iMax,jMin,jMax,
0520 I uIceC, vIceC,
0521 O gUmom, gVmom,
0522 I myTime, myIter, myThid )
4a08d54d3a Mart*0523 DO j=jMin,jMax
0524 DO i=jMin,jMax
0525 rhsU(i,j,bi,bj) = rhsU(i,j,bi,bj) + gUmom(i,j)
0526 rhsV(i,j,bi,bj) = rhsV(i,j,bi,bj) + gVmom(i,j)
8ce59fd29e Mart*0527 ENDDO
0528 ENDDO
0529 ENDIF
0530 #endif /* SEAICE_ALLOW_MOM_ADVECTION */
b7c2bb63b7 Mart*0531 ENDDO
0532 ENDDO
d681af2948 Mart*0533
0534
0535
0536
0537
0538
0539
0540
cbf501ab81 Jean*0541 CALL SEAICE_LSR_CALC_COEFFS(
df1dac8b7b Mart*0542 I etaPlusZeta, zetaMinusEta, etaZ, zetaZ, dragSym,
d681af2948 Mart*0543 O AU, BU, CU, AV, BV, CV, uRt1, uRt2, vRt1, vRt2,
0544 I iMin, iMax, jMin, jMax, myTime, myIter, myThid )
589eca759f Mart*0545
5bb179ddc2 Mart*0546 #ifdef SEAICE_ALLOW_SIDEDRAG
0547 IF ( SEAICEsideDrag .NE. 0. _d 0 ) THEN
0548 DO bj=myByLo(myThid),myByHi(myThid)
0549 DO bi=myBxLo(myThid),myBxHi(myThid)
0550 DO j=jMin,jMax
0551 DO i=iMin,iMax
0552 BU(i,j,bi,bj) = BU(i,j,bi,bj)
0553 & + seaiceMaskU(i,j,bi,bj) * sideDragU(i,j,bi,bj)
0554 BV(i,j,bi,bj) = BV(i,j,bi,bj)
0555 & + seaiceMaskV(i,j,bi,bj) * sideDragV(i,j,bi,bj)
0556 ENDDO
0557 ENDDO
0558 ENDDO
0559 ENDDO
0560 ENDIF
0561 #endif /* SEAICE_ALLOW_SIDEDRAG */
0562
fec75090d3 Jean*0563 #ifndef OBCS_UVICE_OLD
b7c2bb63b7 Mart*0564
0565 DO bj=myByLo(myThid),myByHi(myThid)
0566 DO bi=myBxLo(myThid),myBxHi(myThid)
4a08d54d3a Mart*0567 DO j=jMin,jMax
0568 DO i=iMin,iMax
b7c2bb63b7 Mart*0569 IF ( maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj) .EQ. 0. ) THEN
4a08d54d3a Mart*0570 AU(i,j,bi,bj) = ZERO
0571 BU(i,j,bi,bj) = ONE
0572 CU(i,j,bi,bj) = ZERO
0573 uRt1(i,j,bi,bj) = ZERO
0574 uRt2(i,j,bi,bj) = ZERO
0575 rhsU(i,j,bi,bj) = uIce(i,j,bi,bj)
b7c2bb63b7 Mart*0576 ENDIF
fec75090d3 Jean*0577 IF ( maskInC(i,j,bi,bj)*maskInC(i,j-1,bi,bj) .EQ. 0. ) THEN
4a08d54d3a Mart*0578 AV(i,j,bi,bj) = ZERO
0579 BV(i,j,bi,bj) = ONE
0580 CV(i,j,bi,bj) = ZERO
0581 vRt1(i,j,bi,bj) = ZERO
0582 vRt2(i,j,bi,bj) = ZERO
0583 rhsV(i,j,bi,bj) = vIce(i,j,bi,bj)
fec75090d3 Jean*0584 ENDIF
0585 ENDDO
0586 ENDDO
589eca759f Mart*0587 ENDDO
0588 ENDDO
b7c2bb63b7 Mart*0589 #endif /* OBCS_UVICE_OLD */
589eca759f Mart*0590
7b7a25aa58 Jean*0591
0592
589eca759f Mart*0593 #ifdef ALLOW_DEBUG
86fa82a64d Jean*0594 IF ( debugLevel .GE. debLevD ) THEN
0595 WRITE(msgBuf,'(A,I3,A)')
3e8bf7e743 Mart*0596 & 'Uice pre iter (SEAICE_LSR', MOD(ipass,1000), ')'
7b7a25aa58 Jean*0597 CALL DEBUG_STATS_RL( 1, UICE, msgBuf, myThid )
0598 WRITE(msgBuf,'(A,I3,A)')
3e8bf7e743 Mart*0599 & 'Vice pre iter (SEAICE_LSR', MOD(ipass,1000), ')'
86fa82a64d Jean*0600 CALL DEBUG_STATS_RL( 1, VICE, msgBuf, myThid )
589eca759f Mart*0601 ENDIF
0602 #endif /* ALLOW_DEBUG */
0603
efcb5efb58 Jean*0604
143d9ce879 Dami*0605 IF ( printResidual .OR. LSR_mixIniGuess.GE.1
0606 & .OR. SEAICEuseLSRflex )
0607 & CALL SEAICE_RESIDUAL(
d681af2948 Mart*0608 I rhsU, rhsV, uRt1, uRt2, vRt1, vRt2,
0609 I AU, BU, CU, AV, BV, CV, uIce, vIce,
0610 O residUini, residVini, uTmp, vTmp,
0611 I printResidual, myIter, myThid )
0612
efcb5efb58 Jean*0613 #ifdef SEAICE_ALLOW_FREEDRIFT
d681af2948 Mart*0614 IF ( printResidual .OR. LSR_mixIniGuess.GE.1 ) THEN
0615 IF ( LSR_mixIniGuess.GE.1 ) THEN
0616 DO bj=myByLo(myThid),myByHi(myThid)
0617 DO bi=myBxLo(myThid),myBxHi(myThid)
0618 DO j=jMin,jMax
0619 DO i=iMin,iMax
0620 uIce_fd(i,j,bi,bj) = FORCEX(i,j,bi,bj)
0621 & / ( 1. _d 0 - seaiceMaskU(i,j,bi,bj)
e232688b80 Mart*0622 & + seaiceMassU(i,j,bi,bj)*recip_deltaT
df1dac8b7b Mart*0623 & + HALF*( dragSym(i,j,bi,bj) + dragSym(i-1,j,bi,bj) )
0c560a7486 Mart*0624 & * MAX(areaW(i,j,bi,bj),areaMinLoc)
d681af2948 Mart*0625 & )
0626 vIce_fd(i,j,bi,bj) = FORCEY(i,j,bi,bj)
0627 & / ( 1. _d 0 - seaiceMaskV(i,j,bi,bj)
e232688b80 Mart*0628 & + seaiceMassV(i,j,bi,bj)*recip_deltaT
df1dac8b7b Mart*0629 & + HALF*( dragSym(i,j,bi,bj) + dragSym(i,j-1,bi,bj) )
0c560a7486 Mart*0630 & * MAX(areaS(i,j,bi,bj),areaMinLoc)
d681af2948 Mart*0631 & )
0632 ENDDO
0633 ENDDO
0634 ENDDO
0635 ENDDO
0636 CALL EXCH_UV_XY_RL( uIce_fd, vIce_fd, .TRUE., myThid )
0637 ENDIF
0638 IF ( LSR_mixIniGuess.GE.0 )
efcb5efb58 Jean*0639 & CALL SEAICE_RESIDUAL(
0640 I rhsU, rhsV, uRt1, uRt2, vRt1, vRt2,
0641 I AU, BU, CU, AV, BV, CV, uIce_fd, vIce_fd,
0642 O residU_fd, residV_fd, uRes, vRes,
0643 I printResidual, myIter, myThid )
0644 ENDIF
0645 IF ( LSR_mixIniGuess.GE.2 ) THEN
9b4636ff01 Patr*0646 # ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0647
0648
9b4636ff01 Patr*0649 # endif /* ALLOW_AUTODIFF_TAMC */
efcb5efb58 Jean*0650 DO bj=myByLo(myThid),myByHi(myThid)
0651 DO bi=myBxLo(myThid),myBxHi(myThid)
d681af2948 Mart*0652 DO j=jMin,jMax
0653 DO i=iMin,iMax
0654 errIni = uTmp(i,j,bi,bj)*uTmp(i,j,bi,bj)
0655 errFD = uRes(i,j,bi,bj)*uRes(i,j,bi,bj)
0656 IF ( LSR_mixIniGuess.GE.4 ) THEN
0657 errIni = errIni*errIni
0658 errFD = errFD *errFD
0659 ENDIF
0660 errSum = ( errIni + errFD )
0661 & *maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj)
0662 IF ( errSum.GT.0. ) THEN
0663 uIce(i,j,bi,bj) = ( errFD *uIce(i,j,bi,bj)
0664 & + errIni*uIce_fd(i,j,bi,bj)
0665 & )/errSum
0666 ENDIF
0667 errIni = vTmp(i,j,bi,bj)*vTmp(i,j,bi,bj)
0668 errFD = vRes(i,j,bi,bj)*vRes(i,j,bi,bj)
0669 IF ( LSR_mixIniGuess.GE.4 ) THEN
0670 errIni = errIni*errIni
0671 errFD = errFD *errFD
0672 ENDIF
0673 errSum = ( errIni + errFD )
0674 & *maskInC(i,j,bi,bj)*maskInC(i,j-1,bi,bj)
0675 IF ( errSum.GT.0. ) THEN
0676 vIce(i,j,bi,bj) = ( errFD *vIce(i,j,bi,bj)
0677 & + errIni*vIce_fd(i,j,bi,bj)
0678 & )/errSum
0679 ENDIF
efcb5efb58 Jean*0680 ENDDO
d681af2948 Mart*0681 ENDDO
efcb5efb58 Jean*0682 ENDDO
0683 ENDDO
0684 CALL EXCH_UV_XY_RL( uIce, vIce, .TRUE., myThid )
dfc97bd020 Mart*0685 #ifndef ALLOW_AUTODIFF_TAMC
0e2f719d67 Mart*0686
0687
0688
efcb5efb58 Jean*0689 IF ( printResidual ) THEN
0690 CALL SEAICE_RESIDUAL(
0691 I rhsU, rhsV, uRt1, uRt2, vRt1, vRt2,
0692 I AU, BU, CU, AV, BV, CV, uIce, vIce,
0693 O residUmix, residVmix, uTmp, vTmp,
0694 I printResidual, myIter, myThid )
0695 ENDIF
dfc97bd020 Mart*0696 #endif /* ALLOW_AUTODIFF_TAMC */
efcb5efb58 Jean*0697 ENDIF
0698 #endif /* SEAICE_ALLOW_FREEDRIFT */
0699
0700
0701
589eca759f Mart*0702
0703
7b7a25aa58 Jean*0704 doIterate4u = .TRUE.
0705 doIterate4v = .TRUE.
589eca759f Mart*0706
0707
d681af2948 Mart*0708
cbf501ab81 Jean*0709
964a494b4b Mart*0710
0711
0e2f719d67 Mart*0712
964a494b4b Mart*0713
9b4636ff01 Patr*0714
0d75a51072 Mart*0715 ICOUNT1 = linearIterLoc
0716 ICOUNT2 = linearIterLoc
9b4636ff01 Patr*0717
143d9ce879 Dami*0718 #ifdef SEAICE_ALLOW_LSR_FLEX
0719
0720 ENDIF
0721
0722 IF ( SEAICEuseLSRflex ) THEN
0723
0724 IF ( residUini .EQ. 0. _d 0) THEN
0725 doIterate4u = .FALSE.
0726 ICOUNT1 = 0
0727 S1 = residUini
0728 ENDIF
0729 IF ( residVini .EQ. 0. _d 0) THEN
0730 doIterate4v = .FALSE.
0731 ICOUNT2 = 0
0732 S2 = residVini
0733 ENDIF
0734
0735 residIni = SQRT(residUini*residUini+residVini*residVini)
0736
0737
0738
0739 IF ( ipass .EQ. 1 ) residIniNonLin = residIni
0740
0741
0742 doNonLinLoop = doNonLinLoop .AND.( doIterate4u .OR. doIterate4v )
0743
0744 doNonLinLoop = doNonLinLoop .AND. .NOT.
0745 & ( ipass .GT. 2 .AND.
0746 & ( residIni .LT. SEAICEnonLinTol*residIniNonlin ) )
0747
0748
0749 IF ( residIni .EQ. 0. _d 0 ) residIni = 1. _d -20
0750
0751
0752
0753 LSRflexFac = 1. _d 0 / (1. _d 0 + ABS( LOG10(residIni) ) )
0754
0755
0756 LSRflexFac = MIN(LSRflexFac, 0.99 _d 0 )
0757
0758
0759
0760
0761
0762
0763 LSR_ERRORfu = residUini*LSRflexFac
0764 LSR_ERRORfv = residVini*LSRflexFac
0765
0766 IF ( printResidual ) THEN
0767 _BEGIN_MASTER( myThid )
0768 WRITE(standardMessageUnit,'(A,1X,I5,1X,F10.6,3(1X,E12.6))')
0769 & ' SEAICE_LSR: ipass, FLEX_FACTOR,'//
0770 & ' LSR_ERRORfu, LSR_ERRORfv, residIni =',
0771 & ipass, LSRflexFac, LSR_ERRORfu, LSR_ERRORfv, residIni
0772 _END_MASTER( myThid )
0773 ENDIF
0774
0775 ENDIF
0776
0777
0778 IF ( doNonLinLoop ) THEN
0779 #endif /* SEAICE_ALLOW_LSR_FLEX */
0780
0d75a51072 Mart*0781 DO m = 1, linearIterLoc
9b4636ff01 Patr*0782
0e2f719d67 Mart*0783 #ifdef ALLOW_AUTODIFF_TAMC
0784 # ifdef SEAICE_LSR_ADJOINT_ITER
edb6656069 Mart*0785 lnkey = (nlkey-1)*SOLV_MAX_FIXED + MIN(m,SOLV_MAX_FIXED)
0e2f719d67 Mart*0786 # else
edb6656069 Mart*0787 lnkey = nlkey
0e2f719d67 Mart*0788 # endif /* SEAICE_LSR_ADJOINT_ITER */
3e8bf7e743 Mart*0789
0790
edb6656069 Mart*0791
0792
0e2f719d67 Mart*0793 #endif /* ALLOW_AUTODIFF_TAMC */
9b4636ff01 Patr*0794
7b7a25aa58 Jean*0795 IF ( doIterate4u .OR. doIterate4v ) THEN
589eca759f Mart*0796
7b7a25aa58 Jean*0797 IF ( useCubedSphereExchange ) THEN
0798 doIterate4u = .TRUE.
0799 doIterate4v = .TRUE.
0800 ENDIF
589eca759f Mart*0801
9b4636ff01 Patr*0802 # ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0803
3e8bf7e743 Mart*0804
0805
0806
edb6656069 Mart*0807
9b4636ff01 Patr*0808 # endif /* ALLOW_AUTODIFF_TAMC */
0809
2d6335c1f7 Jean*0810 #ifdef SEAICE_GLOBAL_3DIAG_SOLVER
0811 IF ( SEAICEuseMultiTileSolver ) THEN
0812
0813 DO bj=myByLo(myThid),myByHi(myThid)
0814 DO bi=myBxLo(myThid),myBxHi(myThid)
0815 IF ( doIterate4u ) THEN
5fb12930f9 Mart*0816 DO j=jMin,jMax
0817 DO i=iMin,iMax
2d6335c1f7 Jean*0818 uTmp(i,j,bi,bj) = uIce(i,j,bi,bj)
0819 rhsUx(i,j,bi,bj) = ( rhsU(i,j,bi,bj)
0820 & + uRt1(i,j,bi,bj)*uIce(i,j-1,bi,bj)
0821 & + uRt2(i,j,bi,bj)*uIce(i,j+1,bi,bj)
0822 & )*seaiceMaskU(i,j,bi,bj)
0823 ENDDO
0824 ENDDO
0825 ENDIF
0826 IF ( doIterate4v ) THEN
5fb12930f9 Mart*0827 DO j=jMin,jMax
0828 DO i=iMin,iMax
2d6335c1f7 Jean*0829 vTmp(i,j,bi,bj) = vIce(i,j,bi,bj)
0830 rhsVy(i,j,bi,bj) = ( rhsV(i,j,bi,bj)
4a08d54d3a Mart*0831 & + vRt1(i,j,bi,bj)*vIce(i-1,j,bi,bj)
0832 & + vRt2(i,j,bi,bj)*vIce(i+1,j,bi,bj)
2d6335c1f7 Jean*0833 & )*seaiceMaskU(i,j,bi,bj)
0834 ENDDO
0835 ENDDO
0836 ENDIF
0837
0838 ENDDO
0839 ENDDO
0840
0d75a51072 Mart*0841 iSubIter = linearIterLoc*(ipass-1) + m
2d6335c1f7 Jean*0842 CALL SOLVE_UV_TRIDIAGO(
0843 I 1, OLx, doIterate4u, doIterate4v,
0844 I AU, BU, CU, rhsUx,
0845 I AV, BV, CV, rhsVy,
0846 O uIce, vIce, errCode,
0847 I iSubIter, myIter, myThid )
0848
0849 DO bj=myByLo(myThid),myByHi(myThid)
0850 DO bi=myBxLo(myThid),myBxHi(myThid)
0851 IF ( doIterate4u ) THEN
5fb12930f9 Mart*0852 DO j=jMin,jMax
0853 DO i=iMin,iMax
2d6335c1f7 Jean*0854 uIce(i,j,bi,bj) = uTmp(i,j,bi,bj)
0855 & + WFAU*( uIce(i,j,bi,bj)-uTmp(i,j,bi,bj) )
0856 ENDDO
0857 ENDDO
0858 ENDIF
0859 IF ( doIterate4v ) THEN
5fb12930f9 Mart*0860 DO j=jMin,jMax
0861 DO i=iMin,iMax
2d6335c1f7 Jean*0862 vIce(i,j,bi,bj) = vTmp(i,j,bi,bj)
0863 & + WFAV*( vIce(i,j,bi,bj)-vTmp(i,j,bi,bj) )
0864 ENDDO
0865 ENDDO
0866 ENDIF
0867
0868 ENDDO
0869 ENDDO
0870
0871 ELSE
0872 #endif /* SEAICE_GLOBAL_3DIAG_SOLVER */
0873
7b7a25aa58 Jean*0874 DO bj=myByLo(myThid),myByHi(myThid)
0875 DO bi=myBxLo(myThid),myBxHi(myThid)
0876
f4df908d9c Mart*0877 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0878 tkey = bi + (bj-1)*nSx + (lnkey-1)*nSx*nSy
0879
0880
f4df908d9c Mart*0881 #endif /* ALLOW_AUTODIFF_TAMC */
4a08d54d3a Mart*0882
3e8bf7e743 Mart*0883
8df0ccd026 Jean*0884 DO j=1-OLy,sNy+OLy
0885 DO i=1-OLx,sNx+OLx
4a08d54d3a Mart*0886 uTmp(i,j,bi,bj)=uIce(i,j,bi,bj)
0887 vTmp(i,j,bi,bj)=vIce(i,j,bi,bj)
8df0ccd026 Jean*0888 ENDDO
0889 ENDDO
0890
7b7a25aa58 Jean*0891 IF ( doIterate4u ) THEN
0892
5fb12930f9 Mart*0893 CALL SEAICE_LSR_TRIDIAGU(
8df0ccd026 Jean*0894 I AU, BU, CU, uRt1, uRt2, rhsU, uTmp, seaiceMaskU, WFAU,
5fb12930f9 Mart*0895 U uIce,
8d9be9cde5 Mart*0896 I iMin, iMax, jMin, jMax, bi, bj, myTime, myIter, myThid )
3ff8c9e192 Jean*0897 ENDIF
589eca759f Mart*0898
7b7a25aa58 Jean*0899 IF ( doIterate4v ) THEN
0900
5fb12930f9 Mart*0901 CALL SEAICE_LSR_TRIDIAGV(
cbf501ab81 Jean*0902 I AV, BV, CV, vRt1, vRt2, rhsV, vTmp, seaiceMaskV, WFAV,
5fb12930f9 Mart*0903 U vIce,
8d9be9cde5 Mart*0904 I iMin, iMax, jMin, jMax, bi, bj, myTime, myIter, myThid )
7b7a25aa58 Jean*0905 ENDIF
589eca759f Mart*0906
7b7a25aa58 Jean*0907
3ff8c9e192 Jean*0908 ENDDO
589eca759f Mart*0909 ENDDO
0910
2d6335c1f7 Jean*0911 #ifdef SEAICE_GLOBAL_3DIAG_SOLVER
0912 ENDIF
0913 #endif /* SEAICE_GLOBAL_3DIAG_SOLVER */
0914
143d9ce879 Dami*0915 #ifdef SEAICE_ALLOW_LSR_FLEX
0916 IF ( SEAICEuseLSRflex ) THEN
0917 IF ( MOD(m,SOLV_NCHECK) .EQ. 0
0918 & .AND. (doIterate4u.OR.doIterate4v) ) THEN
0919 CALL SEAICE_RESIDUAL(
0920 I rhsU, rhsV, uRt1, uRt2, vRt1, vRt2,
0921 I AU, BU, CU, AV, BV, CV, uIce, vIce,
0922 O S1, S2, uTmp, vTmp,
0923 I printResidual, myIter, myThid )
0924 IF( S1.LT.LSR_ERRORfu ) THEN
0925 ICOUNT1=m
0926 doIterate4u = .FALSE.
0927 ENDIF
0928 IF( S2.LT.LSR_ERRORfv ) THEN
0929 ICOUNT2=m
0930 doIterate4v = .FALSE.
0931 ENDIF
0932 ENDIF
0933 ELSE
0934 #endif
7b7a25aa58 Jean*0935 IF ( doIterate4u.AND.MOD(m,SOLV_NCHECK).EQ.0) THEN
0936 S1=ZERO
0937 DO bj=myByLo(myThid),myByHi(myThid)
0938 DO bi=myBxLo(myThid),myBxHi(myThid)
4a08d54d3a Mart*0939 DO j=1,sNy
0940 DO i=1,sNx
0941 UERR=(uIce(i,j,bi,bj)-uTmp(i,j,bi,bj))
0942 & * seaiceMaskU(i,j,bi,bj)
7b7a25aa58 Jean*0943 S1=MAX(ABS(UERR),S1)
0944 ENDDO
0945 ENDDO
3ff8c9e192 Jean*0946 ENDDO
0947 ENDDO
7b7a25aa58 Jean*0948 _GLOBAL_MAX_RL( S1, myThid )
0949
3e8bf7e743 Mart*0950
589eca759f Mart*0951
7b7a25aa58 Jean*0952 IF(m.GT.1.AND.S1.GT.S1A) WFAU=WFAU2
0953 S1A=S1
0954 IF(S1.LT.LSR_ERROR) THEN
0955 ICOUNT1=m
0956 doIterate4u = .FALSE.
0957 ENDIF
0958 ENDIF
589eca759f Mart*0959
7b7a25aa58 Jean*0960 IF ( doIterate4v.AND.MOD(m,SOLV_NCHECK).EQ.0) THEN
0961 S2=ZERO
0962 DO bj=myByLo(myThid),myByHi(myThid)
0963 DO bi=myBxLo(myThid),myBxHi(myThid)
4a08d54d3a Mart*0964 DO j=1,sNy
0965 DO i=1,sNx
0966 UERR=(vIce(i,j,bi,bj)-vTmp(i,j,bi,bj))
0967 & * seaiceMaskV(i,j,bi,bj)
7b7a25aa58 Jean*0968 S2=MAX(ABS(UERR),S2)
0969 ENDDO
0970 ENDDO
0971 ENDDO
0972 ENDDO
0973 _GLOBAL_MAX_RL( S2, myThid )
0974
0975 IF(m.GT.1.AND.S2.GT.S2A) WFAV=WFAV2
0976 S2A=S2
0977 IF(S2.LT.LSR_ERROR) THEN
0978 ICOUNT2=m
0979 doIterate4v = .FALSE.
0980 ENDIF
0981 ENDIF
143d9ce879 Dami*0982 #ifdef SEAICE_ALLOW_LSR_FLEX
0983 ENDIF
0984 #endif
7b7a25aa58 Jean*0985
0986 CALL EXCH_UV_XY_RL( uIce, vIce,.TRUE.,myThid)
0987
0988
0989 ENDIF
589eca759f Mart*0990 ENDDO
0991
0992
143d9ce879 Dami*0993 #ifdef SEAICE_ALLOW_LSR_FLEX
0994 IF ( SEAICEuseLSRflex ) THEN
0995 residkm2 = residkm1
0996 residkm1 = residIni
0997 ENDIF
0998 #endif /* SEAICE_ALLOW_LSR_FLEX */
0999
efcb5efb58 Jean*1000
1001
dfc97bd020 Mart*1002 #ifndef ALLOW_AUTODIFF_TAMC
0e2f719d67 Mart*1003
1004
1005
efcb5efb58 Jean*1006 IF ( printResidual ) THEN
1007
1008 CALL SEAICE_RESIDUAL(
1009 I rhsU, rhsV, uRt1, uRt2, vRt1, vRt2,
1010 I AU, BU, CU, AV, BV, CV, uIce, vIce,
1011 O residUend, residVend, uTmp, vTmp,
1012 I printResidual, myIter, myThid )
143d9ce879 Dami*1013
86fa82a64d Jean*1014 _BEGIN_MASTER( myThid )
c512e371cc drin*1015 WRITE(standardMessageUnit,'(A,1X,I9,1P3E16.8)')
e116529eb1 Mart*1016 & ' SEAICE_LSR: Residual Initial ipass,Uice,Vice=',
3e8bf7e743 Mart*1017 & ipass, residUini, residVini
efcb5efb58 Jean*1018 #ifdef SEAICE_ALLOW_FREEDRIFT
1019 IF ( LSR_mixIniGuess.GE.0 )
1020 & WRITE(standardMessageUnit,'(A,1P3E16.8)')
1021 & ' SEAICE_LSR: Residual FrDrift U_fd,V_fd=',
1022 & residU_fd, residV_fd
1023 IF ( LSR_mixIniGuess.GE.2 )
1024 & WRITE(standardMessageUnit,'(A,1P3E16.8)')
1025 & ' SEAICE_LSR: Residual Mix-ini Uice,Vice=',
1026 & residUmix, residVmix
1027 #endif /* SEAICE_ALLOW_FREEDRIFT */
c512e371cc drin*1028 WRITE(standardMessageUnit,'(A,I9,A,I9,1P2E16.8)')
3e8bf7e743 Mart*1029 & ' SEAICE_LSR (ipass=',ipass,') iters,dU,Resid=',
efcb5efb58 Jean*1030 & ICOUNT1, S1, residUend
c512e371cc drin*1031 WRITE(standardMessageUnit,'(A,I9,A,I9,1P2E16.8)')
3e8bf7e743 Mart*1032 & ' SEAICE_LSR (ipass=',ipass,') iters,dV,Resid=',
efcb5efb58 Jean*1033 & ICOUNT2, S2, residVend
86fa82a64d Jean*1034 _END_MASTER( myThid )
7b7a25aa58 Jean*1035 ENDIF
efcb5efb58 Jean*1036 #ifdef ALLOW_DEBUG
7b7a25aa58 Jean*1037 IF ( debugLevel .GE. debLevD ) THEN
1038 WRITE(msgBuf,'(A,I3,A)')
3e8bf7e743 Mart*1039 & 'Uice post iter (SEAICE_LSR', MOD(ipass,1000), ')'
7b7a25aa58 Jean*1040 CALL DEBUG_STATS_RL( 1, UICE, msgBuf, myThid )
1041 WRITE(msgBuf,'(A,I3,A)')
3e8bf7e743 Mart*1042 & 'Vice post iter (SEAICE_LSR', MOD(ipass,1000), ')'
7b7a25aa58 Jean*1043 CALL DEBUG_STATS_RL( 1, VICE, msgBuf, myThid )
589eca759f Mart*1044 ENDIF
1045 #endif /* ALLOW_DEBUG */
bc6ed4e27a Jean*1046 IF ( doIterate4u .OR. doIterate4v ) THEN
1047 WRITE(msgBuf,'(2A,I10,A,I4,A)') '** WARNING ** SEAICE_LSR ',
3e8bf7e743 Mart*1048 & '(it=', myIter, ',', ipass,') did not converge :'
bc6ed4e27a Jean*1049 CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
1050 & SQUEEZE_RIGHT, myThid )
1051 WRITE(msgBuf,'(2(A,I6,0PF6.3,1PE14.6))')
1052 & ' nIt,wFU,dU=', ICOUNT1, WFAU, S1,
1053 & ' ; nIt,wFV,dV=', ICOUNT2, WFAV, S2
1054 CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
1055 & SQUEEZE_RIGHT, myThid )
1056 ENDIF
dfc97bd020 Mart*1057 #endif /* ALLOW_AUTODIFF_TAMC */
589eca759f Mart*1058
1059
1060 DO bj=myByLo(myThid),myByHi(myThid)
1061 DO bi=myBxLo(myThid),myBxHi(myThid)
4a08d54d3a Mart*1062 DO j=1-OLy,sNy+OLy
1063 DO i=1-OLx,sNx+OLx
1064 uIce(i,j,bi,bj)=uIce(i,j,bi,bj)* seaiceMaskU(i,j,bi,bj)
1065 vIce(i,j,bi,bj)=vIce(i,j,bi,bj)* seaiceMaskV(i,j,bi,bj)
3ff8c9e192 Jean*1066 ENDDO
1067 ENDDO
589eca759f Mart*1068 ENDDO
1069 ENDDO
772590b63c Mart*1070
589eca759f Mart*1071
992ae64189 Mart*1072 #ifdef SEAICE_ALLOW_CHECK_LSR_CONVERGENCE
86fa82a64d Jean*1073 IF ( debugLevel .GE. debLevC ) THEN
992ae64189 Mart*1074 resnorm = 0. _d 0
1075 EKnorm = 0. _d 0
1076 counter = 0. _d 0
1077 DO bj=myByLo(myThid),myByHi(myThid)
1078 DO bi=myBxLo(myThid),myBxHi(myThid)
1079
1080 DO j=1,sNy
1081 DO i=1,sNx
1082 resnorm = resnorm + 0.5 _d 0 *
1083 & ( ( (uIceNm1(i,j,bi,bj)+uIceNm1(i+1,j,bi,bj))
1084 & - (uIce(i,j,bi,bj)+uIce(i+1,j,bi,bj)) )**2
1085 & + ( (vIceNm1(i,j,bi,bj)+vIceNm1(i,j+1,bi,bj))
1086 & - (vIce(i,j,bi,bj)+vIce(i,j+1,bi,bj)) )**2 )
8377b8ee87 Mart*1087 IF ( area(i,j,bi,bj) .GT. 0.5 _d 0 ) THEN
992ae64189 Mart*1088 EKnorm = EKnorm + 0.5 _d 0 * heff(i,j,bi,bj) *
1089 & ( ( (uIce(i,j,bi,bj)+uIce(i+1,j,bi,bj)) )**2
1090 & + ( (vIce(i,j,bi,bj)+vIce(i,j+1,bi,bj)) )**2 )
1091 counter = counter + 1. _d 0
1092 ENDIF
1093 ENDDO
1094 ENDDO
1095 ENDDO
1096 ENDDO
1097 _GLOBAL_SUM_RL( resnorm, myThid )
1098 _GLOBAL_SUM_RL( EKnorm, myThid )
1099 _GLOBAL_SUM_RL( counter, myThid )
8377b8ee87 Mart*1100 IF ( counter .GT. 0. _d 0 ) EKnorm = EKnorm/counter
86fa82a64d Jean*1101 _BEGIN_MASTER( myThid )
1102 WRITE(standardMessageUnit,'(A,I7,1X,2E22.14)')
1103 & 'S/R SEAICE_LSR: IPSEUDO, RESNORM, EKNORM = ',
8377b8ee87 Mart*1104 & ipass, SQRT(resnorm), EKnorm
86fa82a64d Jean*1105 _END_MASTER( myThid )
992ae64189 Mart*1106 ENDIF
1107 #endif /* SEAICE_ALLOW_CHECK_LSR_CONVERGENCE */
1108
143d9ce879 Dami*1109
992ae64189 Mart*1110 ENDIF
3e8bf7e743 Mart*1111
992ae64189 Mart*1112 ENDDO
1113
143d9ce879 Dami*1114 #ifdef ALLOW_DIAGNOSTICS
1115 IF ( DIAGNOSTICS_IS_ON('SIlsrRe ',myThid) ) THEN
1116
1117
1118
1119 CALL EXCH_UV_XY_RL( uTmp, vTmp,.TRUE.,myThid)
1120 DO bj = myByLo(myThid), myByHi(myThid)
1121 DO bi = myBxLo(myThid), myBxHi(myThid)
1122 DO j=1,sNy
1123 DO i=1,sNx
1124 RTmp(i,j) = SQRT( HALF * (
1125 & uTmp(i, j,bi,bj) * uTmp(i, j,bi,bj)
1126 & + uTmp(i+1,j,bi,bj) * uTmp(i+1,j,bi,bj)
1127 & + vTmp(i,j, bi,bj) * vTmp(i,j, bi,bj)
1128 & + vTmp(i,j+1,bi,bj) * vTmp(i,j+1,bi,bj)
1129 & ) )
1130 ENDDO
1131 ENDDO
1132 CALL DIAGNOSTICS_FILL(RTmp,'SIlsrRe ',0,1,3,bi,bj,myThid)
1133 ENDDO
1134 ENDDO
1135 ENDIF
1136 #endif /* ALLOW_DIAGNOSTICS */
1137
992ae64189 Mart*1138 IF ( useHB87StressCoupling ) THEN
dfc97bd020 Mart*1139 # ifdef ALLOW_AUTODIFF_TAMC
1140
1141
cbf501ab81 Jean*1142
dfc97bd020 Mart*1143
1144 # endif /* ALLOW_AUTODIFF_TAMC */
8b339058e0 Mart*1145
35fda33b05 Jean*1146
8b339058e0 Mart*1147
1148 CALL SEAICE_CALC_STRAINRATES(
1149 I uIce, vIce,
1150 O e11, e22, e12,
2e75568507 Mart*1151 I 3, myTime, myIter, myThid )
8b339058e0 Mart*1152
925df4f4b9 Mart*1153
8b339058e0 Mart*1154 DO bj=myByLo(myThid),myByHi(myThid)
1155 DO bi=myBxLo(myThid),myBxHi(myThid)
cbf501ab81 Jean*1156 CALL SEAICE_CALC_STRESSDIV(
925df4f4b9 Mart*1157 I e11, e22, e12, press, zeta, eta, etaZ,
1158 O stressDivergenceX, stressDivergenceY,
1159 I bi, bj, myTime, myIter, myThid )
8b339058e0 Mart*1160 ENDDO
964a494b4b Mart*1161 ENDDO
35fda33b05 Jean*1162
8b339058e0 Mart*1163 ENDIF
1164
53092bcb42 Mart*1165 RETURN
1166 END
efcb5efb58 Jean*1167
1168
1169
1170
1171
1172 SUBROUTINE SEAICE_RESIDUAL(
1173 I rhsU, rhsV, uRt1, uRt2, vRt1, vRt2,
1174 I AU, BU, CU, AV, BV, CV, uFld, vFld,
1175 O residU, residV, uRes, vRes,
1176 I calcMeanResid, myIter, myThid )
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186 IMPLICIT NONE
1187
1188
1189 #include "SIZE.h"
1190 #include "EEPARAMS.h"
1191 #include "GRID.h"
03c669d1ab Jean*1192 #include "SEAICE_SIZE.h"
efcb5efb58 Jean*1193 #include "SEAICE.h"
1194
1195
1196
1197
1198
1199
1200
03c669d1ab Jean*1201 _RL rhsU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1202 _RL rhsV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
efcb5efb58 Jean*1203
03c669d1ab Jean*1204 _RL uRt1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1205 _RL uRt2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
efcb5efb58 Jean*1206
03c669d1ab Jean*1207 _RL vRt1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1208 _RL vRt2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
efcb5efb58 Jean*1209
03c669d1ab Jean*1210 _RL AU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1211 _RL BU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1212 _RL CU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1213 _RL AV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1214 _RL BV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1215 _RL CV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
efcb5efb58 Jean*1216
03c669d1ab Jean*1217 _RL uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1218 _RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
efcb5efb58 Jean*1219
1220 _RL residU, residV
03c669d1ab Jean*1221 _RL uRes (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1222 _RL vRes (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
efcb5efb58 Jean*1223 LOGICAL calcMeanResid
1224 INTEGER myIter
1225 INTEGER myThid
1226
1227
1228
1229
1230
1231
1232 INTEGER i, j, bi, bj
1233
1234
1235
1236 residU = 0.
1237 residV = 0.
1238
1239 DO bj=myByLo(myThid),myByHi(myThid)
1240 DO bi=myBxLo(myThid),myBxHi(myThid)
1241 DO j=1,sNy
1242 DO i=1,sNx
1243 uRes(i,j,bi,bj) = rhsU(i,j,bi,bj)
1244 & + uRt1(i,j,bi,bj)*uFld(i,j-1,bi,bj)
1245 & + uRt2(i,j,bi,bj)*uFld(i,j+1,bi,bj)
1246 & - ( AU(i,j,bi,bj)*uFld(i-1,j,bi,bj)
1247 & + BU(i,j,bi,bj)*uFld( i ,j,bi,bj)
1248 & + CU(i,j,bi,bj)*uFld(i+1,j,bi,bj)
1249 & )
1250 vRes(i,j,bi,bj) = rhsV(i,j,bi,bj)
1251 & + vRt1(i,j,bi,bj)*vFld(i-1,j,bi,bj)
1252 & + vRt2(i,j,bi,bj)*vFld(i+1,j,bi,bj)
1253 & - ( AV(i,j,bi,bj)*vFld(i,j-1,bi,bj)
1254 & + BV(i,j,bi,bj)*vFld(i, j ,bi,bj)
1255 & + CV(i,j,bi,bj)*vFld(i,j+1,bi,bj)
1256 & )
1257 ENDDO
1258 ENDDO
56d0d0c8e2 Mart*1259 IF ( calcMeanResid ) THEN
1260 DO j=1,sNy
1261 DO i=1,sNx
efcb5efb58 Jean*1262 residU = residU + uRes(i,j,bi,bj)*uRes(i,j,bi,bj)
1263 & *rAw(i,j,bi,bj)*maskInW(i,j,bi,bj)
5adf1d4e76 Mart*1264 #ifndef OBCS_UVICE_OLD
1265 & *maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj)
1266 #endif
efcb5efb58 Jean*1267 residV = residV + vRes(i,j,bi,bj)*vRes(i,j,bi,bj)
1268 & *rAs(i,j,bi,bj)*maskInS(i,j,bi,bj)
5adf1d4e76 Mart*1269 #ifndef OBCS_UVICE_OLD
1270 & *maskInC(i,j,bi,bj)*maskInC(i,j-1,bi,bj)
1271 #endif
56d0d0c8e2 Mart*1272 ENDDO
efcb5efb58 Jean*1273 ENDDO
56d0d0c8e2 Mart*1274 ENDIF
efcb5efb58 Jean*1275 ENDDO
1276 ENDDO
1277 IF ( calcMeanResid ) THEN
56d0d0c8e2 Mart*1278 _GLOBAL_SUM_RL( residU, myThid )
1279 _GLOBAL_SUM_RL( residV, myThid )
1280
1281 IF ( residU.GT.0. ) residU = SQRT(residU/globalArea)
1282 IF ( residV.GT.0. ) residV = SQRT(residV/globalArea)
efcb5efb58 Jean*1283 ENDIF
1284
1285
1286
1287 RETURN
1288 END
b7c2bb63b7 Mart*1289
1290
5fb12930f9 Mart*1291
b7c2bb63b7 Mart*1292
cbf501ab81 Jean*1293 SUBROUTINE SEAICE_LSR_CALC_COEFFS(
df1dac8b7b Mart*1294 I etaPlusZeta, zetaMinusEta, etaZloc, zetaZloc, dragSym,
b7c2bb63b7 Mart*1295 O AU, BU, CU, AV, BV, CV, uRt1, uRt2, vRt1, vRt2,
d585a5aee9 Mart*1296 I iMin, iMax, jMin, jMax, myTime, myIter, myThid )
b7c2bb63b7 Mart*1297
1298
1299
5fb12930f9 Mart*1300
b7c2bb63b7 Mart*1301
1302
1303
1304
1305
1306 IMPLICIT NONE
1307
1308
1309 #include "SIZE.h"
1310 #include "EEPARAMS.h"
1311 #include "PARAMS.h"
1312 #include "GRID.h"
1313 #include "SEAICE_SIZE.h"
1314 #include "SEAICE_PARAMS.h"
1315 #include "SEAICE.h"
1316
1317
1318
1319
1320
1321
1322 _RL myTime
1323 INTEGER myIter
1324 INTEGER myThid
d585a5aee9 Mart*1325 INTEGER iMin, iMax, jMin, jMax
b7c2bb63b7 Mart*1326
1327
1328 _RL etaPlusZeta (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1329 _RL zetaMinusEta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0fd6b2de1a Mart*1330
1331 _RL etaZloc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1332 _RL zetaZloc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
df1dac8b7b Mart*1333
1334 _RL dragSym (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
b7c2bb63b7 Mart*1335
1336
1337
1338 _RL AU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1339 _RL BU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1340 _RL CU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1341 _RL AV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1342 _RL BV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1343 _RL CV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1344
1345 _RL uRt1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1346 _RL uRt2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1347
1348 _RL vRt1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1349 _RL vRt2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1350
1351
1352
1353
1354
1355
1356 INTEGER i, j, bi, bj
1357 _RL hFacM, hFacP
e232688b80 Mart*1358
1359 _RL bdfAlphaOverDt
0fd6b2de1a Mart*1360
1361 _RL strImpCplFac
b7c2bb63b7 Mart*1362
70e078b38a Mart*1363
1364 _RL areaW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1365 _RL areaS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
b7c2bb63b7 Mart*1366
1367
1368
1369
1370
1371
bfe6f7447e Mart*1372 _RL UXX (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1373 _RL UYY (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1374 _RL UXM (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1375 _RL UYM (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1376 _RL VXX (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1377 _RL VYY (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1378 _RL VXM (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1379 _RL VYM (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
b7c2bb63b7 Mart*1380
1381
e501eee760 Mart*1382
e232688b80 Mart*1383 bdfAlphaOverDt = 1. _d 0
e501eee760 Mart*1384 IF ( SEAICEuseBDF2 ) THEN
1385 IF ( myIter.EQ.nIter0 .AND. SEAICEmomStartBDF.EQ.0 ) THEN
e232688b80 Mart*1386 bdfAlphaOverDt = 1. _d 0
6cbc659de0 Mart*1387 ELSE
e232688b80 Mart*1388 bdfAlphaOverDt = 1.5 _d 0
6cbc659de0 Mart*1389 ENDIF
1390 ENDIF
e232688b80 Mart*1391
1392 bdfAlphaOverDt = bdfAlphaOverDt / SEAICE_deltaTdyn
1393
0fd6b2de1a Mart*1394 strImpCplFac = 0. _d 0
1395 IF ( SEAICEuseStrImpCpl ) strImpCplFac = 1. _d 0
1396
b7c2bb63b7 Mart*1397 DO bj=myByLo(myThid),myByHi(myThid)
1398 DO bi=myBxLo(myThid),myBxHi(myThid)
70e078b38a Mart*1399 IF ( SEAICEscaleSurfStress ) THEN
4a08d54d3a Mart*1400 DO j=jMin,jMax
1401 DO i=iMin,iMax
1402 areaW(i,j) = 0.5 _d 0*(AREA(i,j,bi,bj)+AREA(i-1,j,bi,bj))
1403 areaS(i,j) = 0.5 _d 0*(AREA(i,j,bi,bj)+AREA(i,j-1,bi,bj))
70e078b38a Mart*1404 ENDDO
1405 ENDDO
1406 ELSE
4a08d54d3a Mart*1407 DO j=jMin,jMax
1408 DO i=iMin,iMax
1409 areaW(i,j) = 1. _d 0
1410 areaS(i,j) = 1. _d 0
70e078b38a Mart*1411 ENDDO
1412 ENDDO
1413 ENDIF
b7c2bb63b7 Mart*1414
4a08d54d3a Mart*1415
b7c2bb63b7 Mart*1416
4a08d54d3a Mart*1417 DO j=jMin,jMax
1418 DO i=iMin-1,iMax
b7c2bb63b7 Mart*1419
4a08d54d3a Mart*1420 UXX(i,j) = _dyF(i,j,bi,bj) * etaPlusZeta(i,j,bi,bj)
1421 & * _recip_dxF(i,j,bi,bj)
b7c2bb63b7 Mart*1422
4a08d54d3a Mart*1423 UXM(i,j) = _dyF(i,j,bi,bj) * zetaMinusEta(i,j,bi,bj)
1424 & * k1AtC(i,j,bi,bj) * 0.5 _d 0
b7c2bb63b7 Mart*1425 ENDDO
1426 ENDDO
4a08d54d3a Mart*1427 DO j=jMin,jMax+1
1428 DO i=iMin,iMax
b7c2bb63b7 Mart*1429
4a08d54d3a Mart*1430 UYY(i,j) = _dxV(i,j,bi,bj) *
1431 & ( etaZloc(i,j,bi,bj) + strImpCplFac*zetaZloc(i,j,bi,bj) )
1432 & * _recip_dyU(i,j,bi,bj)
b7c2bb63b7 Mart*1433
4a08d54d3a Mart*1434 UYM(i,j) = _dxV(i,j,bi,bj) * etaZloc(i,j,bi,bj)
1435 & * k2AtZ(i,j,bi,bj) * 0.5 _d 0
b7c2bb63b7 Mart*1436 ENDDO
1437 ENDDO
4a08d54d3a Mart*1438 DO j=jMin,jMax
1439 DO i=iMin,iMax+1
b7c2bb63b7 Mart*1440
4a08d54d3a Mart*1441 VXX(i,j) = _dyU(i,j,bi,bj) *
1442 & ( etaZloc(i,j,bi,bj) + strImpCplFac*zetaZloc(i,j,bi,bj) )
1443 & * _recip_dxV(i,j,bi,bj)
b7c2bb63b7 Mart*1444
4a08d54d3a Mart*1445 VXM(i,j) = _dyU(i,j,bi,bj) * etaZloc(i,j,bi,bj)
1446 & * k1AtZ(i,j,bi,bj) * 0.5 _d 0
b7c2bb63b7 Mart*1447 ENDDO
1448 ENDDO
4a08d54d3a Mart*1449 DO j=jMin-1,jMax
1450 DO i=iMin,iMax
b7c2bb63b7 Mart*1451
4a08d54d3a Mart*1452 VYY(i,j) = _dxF(i,j,bi,bj) * etaPlusZeta(i,j,bi,bj)
1453 & * _recip_dyF(i,j,bi,bj)
b7c2bb63b7 Mart*1454
4a08d54d3a Mart*1455 VYM(i,j) = _dxF(i,j,bi,bj) * zetaMinusEta(i,j,bi,bj)
1456 & * k2AtC(i,j,bi,bj) * 0.5 _d 0
b7c2bb63b7 Mart*1457 ENDDO
1458 ENDDO
1459
1460
1461
1462
1463
4a08d54d3a Mart*1464
1465 DO j=jMin,jMax
1466 DO i=iMin,iMax
1467
1468 AU(i,j,bi,bj)= ( - UXX(i-1,j) + UXM(i-1,j) )
1469 & * seaiceMaskU(i,j,bi,bj)
1470
1471 CU(i,j,bi,bj)= ( - UXX(i ,j) - UXM(i ,j) )
1472 & * seaiceMaskU(i,j,bi,bj)
1473
1474 BU(i,j,bi,bj)=(ONE - seaiceMaskU(i,j,bi,bj)) +
1475 & ( UXX(i-1,j) + UXX(i,j) + UYY(i,j+1) + UYY(i,j)
1476 & + UXM(i-1,j) - UXM(i,j) + UYM(i,j+1) - UYM(i,j)
1477 & ) * seaiceMaskU(i,j,bi,bj)
1478
1479 uRt1(i,j,bi,bj)= UYY(i,j ) + UYM(i,j )
1480
1481 uRt2(i,j,bi,bj)= UYY(i,j+1) - UYM(i,j+1)
b7c2bb63b7 Mart*1482 ENDDO
1483 ENDDO
1484
1485
1486
1487
4a08d54d3a Mart*1488 DO j=jMin,jMax
1489 DO i=iMin,iMax
1490 hFacM = seaiceMaskU(i,j-1,bi,bj)
1491 hFacP = seaiceMaskU(i,j+1,bi,bj)
1492
b7c2bb63b7 Mart*1493
1494
4a08d54d3a Mart*1495 BU(i,j,bi,bj)=BU(i,j,bi,bj) + seaiceMaskU(i,j,bi,bj) *
1496 & ( ( 1. _d 0 - hFacM ) * ( UYY(i ,j ) + UYM(i ,j ) )
1497 & + ( 1. _d 0 - hFacP ) * ( UYY(i ,j+1) - UYM(i ,j+1) ) )
1498
1499 uRt1(i,j,bi,bj) = uRt1(i,j,bi,bj) * hFacM
1500 uRt2(i,j,bi,bj) = uRt2(i,j,bi,bj) * hFacP
b7c2bb63b7 Mart*1501 ENDDO
1502 ENDDO
1503
1504
4a08d54d3a Mart*1505 DO j=jMin,jMax
1506 DO i=iMin,iMax
1507 AU(i,j,bi,bj) = AU(i,j,bi,bj) * recip_rAw(i,j,bi,bj)
1508 CU(i,j,bi,bj) = CU(i,j,bi,bj) * recip_rAw(i,j,bi,bj)
cbf501ab81 Jean*1509
1510
e232688b80 Mart*1511
4a08d54d3a Mart*1512 BU(i,j,bi,bj) = BU(i,j,bi,bj) * recip_rAw(i,j,bi,bj)
1513 & + seaiceMaskU(i,j,bi,bj) *
1514 & ( bdfAlphaOverDt*seaiceMassU(i,j,bi,bj)
8377b8ee87 Mart*1515 & + 0.5 _d 0 * ( dragSym(i, j,bi,bj)
4a08d54d3a Mart*1516 & + dragSym(i-1,j,bi,bj) )*areaW(i,j)
df1dac8b7b Mart*1517 & )
4a08d54d3a Mart*1518 uRt1(i,j,bi,bj) = uRt1(i,j,bi,bj) * recip_rAw(i,j,bi,bj)
1519 uRt2(i,j,bi,bj) = uRt2(i,j,bi,bj) * recip_rAw(i,j,bi,bj)
b7c2bb63b7 Mart*1520 ENDDO
1521 ENDDO
1522
1523
1524
1525
1526
4a08d54d3a Mart*1527
1528 DO j=jMin,jMax
1529 DO i=iMin,iMax
1530
1531 AV(i,j,bi,bj)=( - VYY(i,j-1) + VYM(i,j-1)
1532 & ) * seaiceMaskV(i,j,bi,bj)
1533
1534 CV(i,j,bi,bj)=( - VYY(i,j ) - VYM(i,j )
1535 & ) * seaiceMaskV(i,j,bi,bj)
1536
1537 BV(i,j,bi,bj)= (ONE - seaiceMaskV(i,j,bi,bj)) +
1538 & ( VXX(i,j) + VXX(i+1,j) + VYY(i,j) + VYY(i,j-1)
1539 & - VXM(i,j) + VXM(i+1,j) - VYM(i,j) + VYM(i,j-1)
1540 & ) * seaiceMaskV(i,j,bi,bj)
1541
1542 vRt1(i,j,bi,bj) = VXX(i ,j) + VXM(i ,j)
1543
1544 vRt2(i,j,bi,bj) = VXX(i+1,j) - VXM(i+1,j)
b7c2bb63b7 Mart*1545 ENDDO
1546 ENDDO
1547
1548
1549
1550
4a08d54d3a Mart*1551 DO j=jMin,jMax
1552 DO i=iMin,iMax
b7c2bb63b7 Mart*1553 hFacM = seaiceMaskV(i-1,j,bi,bj)
1554 hFacP = seaiceMaskV(i+1,j,bi,bj)
4a08d54d3a Mart*1555
b7c2bb63b7 Mart*1556
1557
4a08d54d3a Mart*1558 BV(i,j,bi,bj)=BV(i,j,bi,bj) + seaiceMaskV(i,j,bi,bj) *
1559 & ( ( 1. _d 0 - hFacM ) * ( VXX(i ,j) + VXM(i ,j) )
1560 & + ( 1. _d 0 - hFacP ) * ( VXX(i+1,j) - VXM(i+1,j) ) )
1561
1562 vRt1(i,j,bi,bj) = vRt1(i,j,bi,bj) * hFacM
1563 vRt2(i,j,bi,bj) = vRt2(i,j,bi,bj) * hFacP
b7c2bb63b7 Mart*1564 ENDDO
1565 ENDDO
1566
1567
4a08d54d3a Mart*1568 DO j=jMin,jMax
1569 DO i=iMin,iMax
1570 AV(i,j,bi,bj) = AV(i,j,bi,bj) * recip_rAs(i,j,bi,bj)
1571 CV(i,j,bi,bj) = CV(i,j,bi,bj) * recip_rAs(i,j,bi,bj)
e232688b80 Mart*1572
cbf501ab81 Jean*1573
e232688b80 Mart*1574
4a08d54d3a Mart*1575 BV(i,j,bi,bj) = BV(i,j,bi,bj) * recip_rAs(i,j,bi,bj)
1576 & + seaiceMaskV(i,j,bi,bj) *
1577 & ( bdfAlphaOverDt*seaiceMassV(i,j,bi,bj)
1578 & + 0.5 _d 0 * ( dragSym(i,j, bi,bj)
1579 & + dragSym(i,j-1,bi,bj) )*areaS(i,j)
df1dac8b7b Mart*1580 & )
4a08d54d3a Mart*1581 vRt1(i,j,bi,bj) = vRt1(i,j,bi,bj) * recip_rAs(i,j,bi,bj)
1582 vRt2(i,j,bi,bj) = vRt2(i,j,bi,bj) * recip_rAs(i,j,bi,bj)
b7c2bb63b7 Mart*1583 ENDDO
1584 ENDDO
1585
70e078b38a Mart*1586 IF ( ( useCubedSphereExchange .AND.
1587 & (SEAICE_OLx.GT.0 .OR. SEAICE_OLy.GT.0) ) .OR.
1588 & SEAICEscaleSurfStress ) THEN
05180ad1ef Mart*1589
1590
1591
1592
70e078b38a Mart*1593
1594
4a08d54d3a Mart*1595 DO j=jMin,jMax
1596 DO i=iMin,iMax
1597 IF (BU(i,j,bi,bj).EQ.0 _d 0) BU(i,j,bi,bj) = 1. _d 0
1598 IF (BV(i,j,bi,bj).EQ.0 _d 0) BV(i,j,bi,bj) = 1. _d 0
05180ad1ef Mart*1599 ENDDO
1600 ENDDO
1601 ENDIF
b7c2bb63b7 Mart*1602
1603 ENDDO
1604 ENDDO
1605
5fb12930f9 Mart*1606 RETURN
1607 END
1608
1609
bfe6f7447e Mart*1610
1611
1612
1613
1614 SUBROUTINE SEAICE_LSR_RHSU(
cbf501ab81 Jean*1615 I zetaMinusEta, etaPlusZeta, etaZloc, zetaZloc, pressLoc,
1616 I uIceC, vIceC,
bfe6f7447e Mart*1617 U rhsU,
1618 I iMin, iMax, jMin, jMax, bi, bj, myThid )
1619
1620
1621
1622
1623
1624
1625
1626
1627
1628 IMPLICIT NONE
1629 #include "SIZE.h"
0fd6b2de1a Mart*1630 #include "EEPARAMS.h"
bfe6f7447e Mart*1631 #include "GRID.h"
1632 #include "SEAICE_SIZE.h"
0fd6b2de1a Mart*1633 #include "SEAICE_PARAMS.h"
bfe6f7447e Mart*1634 #include "SEAICE.h"
1635
1636
1637
1638
1639 INTEGER myThid
1640 INTEGER iMin, iMax, jMin, jMax, bi, bj
1641
1642 _RL zetaMinusEta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1643 _RL etaPlusZeta (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1644 _RL etaZloc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0fd6b2de1a Mart*1645 _RL zetaZloc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
bfe6f7447e Mart*1646 _RL pressloc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0fd6b2de1a Mart*1647 _RL uIceC (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
bfe6f7447e Mart*1648 _RL vIceC (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1649 _RL rhsU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1650
1651
1652
1653
4a08d54d3a Mart*1654 INTEGER i,j
bfe6f7447e Mart*1655 _RL hFacM
1656 _RL sig11(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1657 _RL sig12(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1658
1659
4a08d54d3a Mart*1660 DO j=1-OLy,sNy+OLy
1661 DO i=1-OLx,sNx+OLx
1662 sig11(i,j) = 0. _d 0
1663 sig12(i,j) = 0. _d 0
bfe6f7447e Mart*1664 ENDDO
1665 ENDDO
1666
4a08d54d3a Mart*1667 DO j=jMin,jMax
1668 DO i=iMin-1,iMax
1669 sig11(i,j) = zetaMinusEta(i,j,bi,bj)
1670 & * ( vIceC(i,j+1,bi,bj) - vIceC(i,j,bi,bj) )
1671 & * _recip_dyF(i,j,bi,bj)
1672 & + etaPlusZeta(i,j,bi,bj) * k2AtC(i,j,bi,bj)
1673 & * 0.5 _d 0 * ( vIceC(i,j+1,bi,bj) + vIceC(i,j,bi,bj) )
1674 & - 0.5 _d 0 * pressLoc(i,j,bi,bj)
bfe6f7447e Mart*1675 ENDDO
1676 ENDDO
1677
4a08d54d3a Mart*1678 DO j=jMin,jMax+1
1679 DO i=iMin,iMax
1680 hFacM = seaiceMaskV(i,j,bi,bj) - seaiceMaskV(i-1,j,bi,bj)
1681 sig12(i,j) = etaZloc(i,j,bi,bj) * (
1682 & ( vIceC(i,j,bi,bj) - vIceC(i-1,j,bi,bj) )
1683 & * _recip_dxV(i,j,bi,bj)
1684 & - k1AtZ(i,j,bi,bj)
1685 & * 0.5 _d 0 * ( vIceC(i,j,bi,bj) + vIceC(i-1,j,bi,bj) )
bfe6f7447e Mart*1686 & )
1687
4a08d54d3a Mart*1688 & *HEFFM(i ,j ,bi,bj)*HEFFM(i-1,j ,bi,bj)
1689 & *HEFFM(i ,j-1,bi,bj)*HEFFM(i-1,j-1,bi,bj)
bfe6f7447e Mart*1690
1691
1692
4a08d54d3a Mart*1693 & + etaZloc(i,j,bi,bj) * _recip_dxV(i,j,bi,bj)
1694 & * ( vIceC(i,j,bi,bj) + vIceC(i-1,j,bi,bj) )
bfe6f7447e Mart*1695 & * hFacM * 2. _d 0
1696 ENDDO
1697 ENDDO
1698
0fd6b2de1a Mart*1699 IF ( SEAICEuseStrImpCpl ) THEN
cbf501ab81 Jean*1700
0fd6b2de1a Mart*1701
1702
4a08d54d3a Mart*1703 DO j=jMin,jMax+1
1704 DO i=iMin,iMax
1705 hFacM = seaiceMaskV(i,j,bi,bj) - seaiceMaskV(i-1,j,bi,bj)
1706 sig12(i,j) = sig12(i,j) - zetaZloc(i,j,bi,bj) * (
1707 & ( uIceC(i,j,bi,bj) - uIceC(i,j-1,bi,bj) )
1708 & * _recip_dyU(i,j,bi,bj)
0fd6b2de1a Mart*1709 & )
1710
4a08d54d3a Mart*1711 & *HEFFM(i ,j ,bi,bj)*HEFFM(i-1,j ,bi,bj)
1712 & *HEFFM(i ,j-1,bi,bj)*HEFFM(i-1,j-1,bi,bj)
0fd6b2de1a Mart*1713
1714
1715
4a08d54d3a Mart*1716 & - zetaZloc(i,j,bi,bj) * _recip_dyU(i,j,bi,bj)
1717 & * ( uIceC(i,j,bi,bj) + uIceC(i,j-1,bi,bj) )
0fd6b2de1a Mart*1718 & * hFacM * 2. _d 0
1719 ENDDO
1720 ENDDO
1721 ENDIF
cbf501ab81 Jean*1722
4a08d54d3a Mart*1723 DO j=jMin,jMax
1724 DO i=iMin,iMax
bfe6f7447e Mart*1725
4a08d54d3a Mart*1726 rhsU(i,j,bi,bj) = rhsU(i,j,bi,bj)
1727 & + recip_rAw(i,j,bi,bj) * seaiceMaskU(i,j,bi,bj) *
1728 & ( _dyF(i ,j ,bi,bj)*sig11(i ,j )
1729 & - _dyF(i-1,j ,bi,bj)*sig11(i-1,j )
1730 & + _dxV(i ,j+1,bi,bj)*sig12(i ,j+1)
1731 & - _dxV(i ,j ,bi,bj)*sig12(i ,j ) )
bfe6f7447e Mart*1732 ENDDO
1733 ENDDO
1734
1735 RETURN
1736 END
1737
1738
1739
1740
1741
1742
1743 SUBROUTINE SEAICE_LSR_RHSV(
0fd6b2de1a Mart*1744 I zetaMinusEta, etaPlusZeta, etaZloc, zetaZloc, pressLoc,
1745 I uIceC, vIceC,
bfe6f7447e Mart*1746 U rhsV,
1747 I iMin, iMax, jMin, jMax, bi, bj, myThid )
1748
1749
1750
1751
1752
1753
1754
1755
1756
1757 IMPLICIT NONE
1758 #include "SIZE.h"
0fd6b2de1a Mart*1759 #include "EEPARAMS.h"
bfe6f7447e Mart*1760 #include "GRID.h"
1761 #include "SEAICE_SIZE.h"
0fd6b2de1a Mart*1762 #include "SEAICE_PARAMS.h"
bfe6f7447e Mart*1763 #include "SEAICE.h"
1764
1765
1766
1767
1768 INTEGER myThid
1769 INTEGER iMin, iMax, jMin, jMax, bi, bj
1770
1771 _RL zetaMinusEta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1772 _RL etaPlusZeta (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1773 _RL etaZloc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0fd6b2de1a Mart*1774 _RL zetaZloc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
bfe6f7447e Mart*1775 _RL pressloc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1776 _RL uIceC (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0fd6b2de1a Mart*1777 _RL vIceC (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
bfe6f7447e Mart*1778 _RL rhsV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1779
1780
1781
1782
4a08d54d3a Mart*1783 INTEGER i,j
bfe6f7447e Mart*1784 _RL hFacM
1785 _RL sig22(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1786 _RL sig12(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1787
1788
4a08d54d3a Mart*1789 DO j=1-OLy,sNy+OLy
1790 DO i=1-OLx,sNx+OLx
1791 sig22(i,j) = 0. _d 0
1792 sig12(i,j) = 0. _d 0
bfe6f7447e Mart*1793 ENDDO
1794 ENDDO
1795
1796
4a08d54d3a Mart*1797 DO j=jMin-1,jMax
1798 DO i=iMin,iMax
1799 sig22(i,j) = zetaMinusEta(i,j,bi,bj)
1800 & * ( uIceC(i+1,j,bi,bj) - uIceC(i,j,bi,bj) )
1801 & * _recip_dxF(i,j,bi,bj)
1802 & + etaPlusZeta(i,j,bi,bj) * k1AtC(i,j,bi,bj)
1803 & * 0.5 _d 0 * ( uIceC(i+1,j,bi,bj) + uIceC(i,j,bi,bj) )
1804 & - 0.5 _d 0 * pressLoc(i,j,bi,bj)
bfe6f7447e Mart*1805 ENDDO
1806 ENDDO
1807
4a08d54d3a Mart*1808 DO j=jMin,jMax
1809 DO i=iMin,iMax+1
bfe6f7447e Mart*1810 hFacM = seaiceMaskU(i,j,bi,bj) - seaiceMaskU(i,j-1,bi,bj)
4a08d54d3a Mart*1811 sig12(i,j) = etaZloc(i,j,bi,bj) * (
1812 & ( uIceC(i,j,bi,bj) - uIceC(i,j-1,bi,bj) )
1813 & * _recip_dyU(i,j,bi,bj)
1814 & - k2AtZ(i,j,bi,bj)
1815 & * 0.5 _d 0 * ( uIceC(i,j,bi,bj) + uIceC(i,j-1,bi,bj) )
bfe6f7447e Mart*1816 & )
1817
4a08d54d3a Mart*1818 & *HEFFM(i ,j ,bi,bj)*HEFFM(i-1,j ,bi,bj)
1819 & *HEFFM(i ,j-1,bi,bj)*HEFFM(i-1,j-1,bi,bj)
bfe6f7447e Mart*1820
1821
1822
4a08d54d3a Mart*1823 & + etaZloc(i,j,bi,bj) * _recip_dyU(i,j,bi,bj)
1824 & * ( uIceC(i,j,bi,bj) + uIceC(i,j-1,bi,bj) )
bfe6f7447e Mart*1825 & * hFacM * 2. _d 0
1826 ENDDO
1827 ENDDO
1828
0fd6b2de1a Mart*1829 IF ( SEAICEuseStrImpCpl ) THEN
cbf501ab81 Jean*1830
0fd6b2de1a Mart*1831
1832
4a08d54d3a Mart*1833 DO j=jMin,jMax
1834 DO i=iMin,iMax+1
0fd6b2de1a Mart*1835 hFacM = seaiceMaskU(i,j,bi,bj) - seaiceMaskU(i,j-1,bi,bj)
4a08d54d3a Mart*1836 sig12(i,j) = sig12(i,j) - zetaZloc(i,j,bi,bj) * (
1837 & ( vIceC(i,j,bi,bj) - vIceC(i-1,j,bi,bj) )
1838 & * _recip_dxV(i,j,bi,bj)
0fd6b2de1a Mart*1839 & )
1840
4a08d54d3a Mart*1841 & *HEFFM(i ,j ,bi,bj)*HEFFM(i-1,j ,bi,bj)
1842 & *HEFFM(i ,j-1,bi,bj)*HEFFM(i-1,j-1,bi,bj)
0fd6b2de1a Mart*1843
1844
1845
4a08d54d3a Mart*1846 & - zetaZloc(i,j,bi,bj) * _recip_dxV(i,j,bi,bj)
1847 & * ( vIceC(i,j,bi,bj) + vIceC(i-1,j,bi,bj) )
0fd6b2de1a Mart*1848 & * hFacM * 2. _d 0
1849 ENDDO
1850 ENDDO
1851 ENDIF
1852
4a08d54d3a Mart*1853 DO j=jMin,jMax
1854 DO i=iMin,iMax
bfe6f7447e Mart*1855
4a08d54d3a Mart*1856 rhsV(i,j,bi,bj) = rhsV(i,j,bi,bj)
1857 & + recip_rAs(i,j,bi,bj) * seaiceMaskV(i,j,bi,bj) *
1858 & ( _dyU(i+1,j ,bi,bj) * sig12(i+1,j )
1859 & - _dyU(i ,j ,bi,bj) * sig12(i ,j )
1860 & + _dxF(i ,j ,bi,bj) * sig22(i ,j )
1861 & - _dxF(i ,j-1,bi,bj) * sig22(i ,j-1) )
bfe6f7447e Mart*1862 ENDDO
1863 ENDDO
1864
1865 RETURN
1866 END
1867
1868
5fb12930f9 Mart*1869
1870
1871
1872
1873 SUBROUTINE SEAICE_LSR_TRIDIAGU(
8df0ccd026 Jean*1874 I AU, BU, CU, uRt1, uRt2, rhsU, uTmp, seaiceMaskU, WFAU,
5fb12930f9 Mart*1875 U uIce,
8d9be9cde5 Mart*1876 I iMin, iMax, jMin, jMax, bi, bj, myTime, myIter, myThid )
5fb12930f9 Mart*1877
1878
1879
8d9be9cde5 Mart*1880
5fb12930f9 Mart*1881
1882
1883
1884
1885
1886 IMPLICIT NONE
1887 #include "SIZE.h"
1888
1889
1890
1891
1892
1893
1894 _RL myTime
1895 INTEGER myIter
1896 INTEGER myThid
1897 INTEGER iMin, iMax, jMin, jMax, bi, bj
1898
1899 _RL AU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1900 _RL BU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1901 _RL CU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1902
1903 _RL uRt1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1904 _RL uRt2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1905
1906 _RL rhsU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
8df0ccd026 Jean*1907
1908 _RL uTmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
5fb12930f9 Mart*1909
1910 _RL uIce(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1911
1912 _RL seaiceMaskU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
1913
1914 _RL WFAU
1915
1916
1917
1918
1919
4a08d54d3a Mart*1920 INTEGER i,j,iM
032b010f0c Mart*1921 INTEGER jMinLoc, jStep
12dace3fd6 Mart*1922 #ifdef SEAICE_LSR_ZEBRA
4a08d54d3a Mart*1923 INTEGER k
032b010f0c Mart*1924 PARAMETER ( jStep = 2 )
1925 #else
1926 PARAMETER ( jStep = 1 )
12dace3fd6 Mart*1927 #endif /* SEAICE_LSR_ZEBRA */
7c0da9bbfa Mart*1928 _RL AA3, bet
5fb12930f9 Mart*1929 _RL URT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1930 _RL CUU(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
cbf501ab81 Jean*1931
dadeed8422 Mart*1932
4a08d54d3a Mart*1933 DO j=1-OLy,sNy+OLy
1934 DO i=1-OLx,sNx+OLx
1935 URT(i,j) = 0. _d 0
1936 CUU(i,j) = 0. _d 0
dadeed8422 Mart*1937 ENDDO
1938 ENDDO
12dace3fd6 Mart*1939 #ifdef SEAICE_LSR_ZEBRA
4a08d54d3a Mart*1940 DO k=0,1
1941 jMinLoc = jMin+k
8d9be9cde5 Mart*1942 #else
032b010f0c Mart*1943 jMinLoc = jMin
12dace3fd6 Mart*1944 #endif /* SEAICE_LSR_ZEBRA */
4a08d54d3a Mart*1945 DO j=jMinLoc,jMax,jStep
1946 DO i=iMin,iMax
5fb12930f9 Mart*1947 AA3 = 0. _d 0
4a08d54d3a Mart*1948 IF (i.EQ.iMin) AA3 = AA3 - AU(i,j,bi,bj)*uIce(i-1,j,bi,bj)
1949 IF (i.EQ.iMax) AA3 = AA3 - CU(i,j,bi,bj)*uIce(i+1,j,bi,bj)
cbf501ab81 Jean*1950
4a08d54d3a Mart*1951 URT(i,j)=rhsU(i,j,bi,bj)
5fb12930f9 Mart*1952 & + AA3
4a08d54d3a Mart*1953 & + uRt1(i,j,bi,bj)*uIce(i,j-1,bi,bj)
1954 & + uRt2(i,j,bi,bj)*uIce(i,j+1,bi,bj)
1955 URT(i,j)=URT(i,j)* seaiceMaskU(i,j,bi,bj)
5fb12930f9 Mart*1956 ENDDO
032b010f0c Mart*1957
4a08d54d3a Mart*1958 DO i=iMin,iMax
1959 CUU(i,j)=CU(i,j,bi,bj)
5fb12930f9 Mart*1960 ENDDO
4a08d54d3a Mart*1961 CUU(iMin,j)=CUU(iMin,j)/BU(iMin,j,bi,bj)
1962 URT(iMin,j)=URT(iMin,j)/BU(iMin,j,bi,bj)
5fb12930f9 Mart*1963 #ifdef SEAICE_VECTORIZE_LSR
1964 ENDDO
1965
032b010f0c Mart*1966
4a08d54d3a Mart*1967 DO i=iMin+1,iMax
1968 iM=i-1
1969 DO j=jMinLoc,jMax,jStep
5fb12930f9 Mart*1970 #else /* do not SEAICE_VECTORIZE_LSR */
032b010f0c Mart*1971
4a08d54d3a Mart*1972 DO i=iMin+1,iMax
1973 iM=i-1
5fb12930f9 Mart*1974 #endif /* SEAICE_VECTORIZE_LSR */
4a08d54d3a Mart*1975 bet = BU(i,j,bi,bj)-AU(i,j,bi,bj)*CUU(iM,j)
1976 CUU(i,j) = CUU(i,j)/bet
1977 URT(i,j) = (URT(i,j)-AU(i,j,bi,bj)*URT(iM,j))/bet
5fb12930f9 Mart*1978 ENDDO
1979 #ifdef SEAICE_VECTORIZE_LSR
1980 ENDDO
032b010f0c Mart*1981
4a08d54d3a Mart*1982 DO i=iMin,iMax-1
8377b8ee87 Mart*1983 iM=sNx-i
4a08d54d3a Mart*1984 DO j=jMinLoc,jMax,jStep
8d9be9cde5 Mart*1985 #else
032b010f0c Mart*1986
4a08d54d3a Mart*1987 DO i=iMin,iMax-1
8377b8ee87 Mart*1988 iM=sNx-i
032b010f0c Mart*1989 #endif /* SEAICE_VECTORIZE_LSR */
4a08d54d3a Mart*1990 URT(iM,j)=URT(iM,j)-CUU(iM,j)*URT(iM+1,j)
5fb12930f9 Mart*1991 ENDDO
032b010f0c Mart*1992 #ifdef SEAICE_VECTORIZE_LSR
1993 ENDDO
1994 #endif /* SEAICE_VECTORIZE_LSR */
1995
1996
1997 #ifdef SEAICE_VECTORIZE_LSR
1998
4a08d54d3a Mart*1999 DO j=jMinLoc,jMax,jStep
032b010f0c Mart*2000 #endif /* SEAICE_VECTORIZE_LSR */
4a08d54d3a Mart*2001 DO i=iMin,iMax
2002 uIce(i,j,bi,bj)=uTmp(i,j,bi,bj)
2003 & +WFAU*(URT(i,j)-uTmp(i,j,bi,bj))
5fb12930f9 Mart*2004 ENDDO
2005 ENDDO
12dace3fd6 Mart*2006 #ifdef SEAICE_LSR_ZEBRA
8d9be9cde5 Mart*2007 ENDDO
12dace3fd6 Mart*2008 #endif /* SEAICE_LSR_ZEBRA */
5fb12930f9 Mart*2009
2010 RETURN
2011 END
2012
2013
2014
2015
2016
2017
2018 SUBROUTINE SEAICE_LSR_TRIDIAGV(
cbf501ab81 Jean*2019 I AV, BV, CV, vRt1, vRt2, rhsV, vTmp, seaiceMaskV, WFAV,
5fb12930f9 Mart*2020 U vIce,
8d9be9cde5 Mart*2021 I iMin, iMax, jMin, jMax, bi, bj, myTime, myIter, myThid )
5fb12930f9 Mart*2022
2023
2024
2025
2026
2027
2028
2029
2030
2031 IMPLICIT NONE
2032 #include "SIZE.h"
2033
2034
2035
2036
2037
2038
2039 _RL myTime
2040 INTEGER myIter
2041 INTEGER myThid
2042 INTEGER iMin, iMax, jMin, jMax, bi, bj
2043
2044 _RL AV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
2045 _RL BV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
2046 _RL CV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
2047
2048 _RL vRt1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
2049 _RL vRt2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
2050
2051 _RL rhsV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
8df0ccd026 Jean*2052
2053 _RL vTmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
5fb12930f9 Mart*2054
2055 _RL vIce(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
2056
2057 _RL seaiceMaskV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
2058
2059 _RL WFAV
2060
2061
2062
2063
2064
4a08d54d3a Mart*2065 INTEGER i,j,jM
032b010f0c Mart*2066 INTEGER iMinLoc, iStep
12dace3fd6 Mart*2067 #ifdef SEAICE_LSR_ZEBRA
4a08d54d3a Mart*2068 INTEGER k
032b010f0c Mart*2069 PARAMETER ( iStep = 2 )
2070 #else
2071 PARAMETER ( iStep = 1 )
12dace3fd6 Mart*2072 #endif /* SEAICE_LSR_ZEBRA */
7c0da9bbfa Mart*2073 _RL AA3, bet
5fb12930f9 Mart*2074 _RL VRT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
2075 _RL CVV(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
2076
dadeed8422 Mart*2077
4a08d54d3a Mart*2078 DO j=1-OLy,sNy+OLy
2079 DO i=1-OLx,sNx+OLx
2080 VRT(i,j) = 0. _d 0
2081 CVV(i,j) = 0. _d 0
dadeed8422 Mart*2082 ENDDO
2083 ENDDO
12dace3fd6 Mart*2084 #ifdef SEAICE_LSR_ZEBRA
4a08d54d3a Mart*2085 DO k=0,1
2086 iMinLoc = iMin+k
032b010f0c Mart*2087 #else
2088 iMinLoc = iMin
12dace3fd6 Mart*2089 #endif /* SEAICE_LSR_ZEBRA */
8d9be9cde5 Mart*2090 #ifdef SEAICE_VECTORIZE_LSR
4a08d54d3a Mart*2091 DO j=jMin,jMax
2092 DO i=iMinLoc,iMax,iStep
8d9be9cde5 Mart*2093 #else
4a08d54d3a Mart*2094 DO i=iMinLoc,iMax,iStep
2095 DO j=jMin,jMax
8d9be9cde5 Mart*2096 #endif /* SEAICE_VECTORIZE_LSR */
5fb12930f9 Mart*2097 AA3 = 0. _d 0
4a08d54d3a Mart*2098 IF (j.EQ.jMin) AA3 = AA3 - AV(i,j,bi,bj)*vIce(i,j-1,bi,bj)
2099 IF (j.EQ.jMax) AA3 = AA3 - CV(i,j,bi,bj)*vIce(i,j+1,bi,bj)
5fb12930f9 Mart*2100
4a08d54d3a Mart*2101 VRT(i,j)=rhsV(i,j,bi,bj)
5fb12930f9 Mart*2102 & + AA3
4a08d54d3a Mart*2103 & + vRt1(i,j,bi,bj)*vIce(i-1,j,bi,bj)
2104 & + vRt2(i,j,bi,bj)*vIce(i+1,j,bi,bj)
2105 VRT(i,j)=VRT(i,j)* seaiceMaskV(i,j,bi,bj)
5fb12930f9 Mart*2106
032b010f0c Mart*2107
4a08d54d3a Mart*2108 CVV(i,j)=CV(i,j,bi,bj)
5fb12930f9 Mart*2109 ENDDO
8d9be9cde5 Mart*2110 #ifdef SEAICE_VECTORIZE_LSR
2111 ENDDO
4a08d54d3a Mart*2112 DO i=iMinLoc,iMax,iStep
8d9be9cde5 Mart*2113 #endif /* SEAICE_VECTORIZE_LSR */
4a08d54d3a Mart*2114 CVV(i,jMin)=CVV(i,jMin)/BV(i,jMin,bi,bj)
2115 VRT(i,jMin)=VRT(i,jMin)/BV(i,jMin,bi,bj)
8d9be9cde5 Mart*2116 #ifdef SEAICE_VECTORIZE_LSR
2117 ENDDO
2118
032b010f0c Mart*2119
4a08d54d3a Mart*2120 DO j=jMin+1,jMax
2121 jM=j-1
2122 DO i=iMinLoc,iMax,iStep
8d9be9cde5 Mart*2123 #else /* do not SEAICE_VECTORIZE_LSR */
032b010f0c Mart*2124
4a08d54d3a Mart*2125 DO j=jMin+1,jMax
2126 jM=j-1
8d9be9cde5 Mart*2127 #endif /* SEAICE_VECTORIZE_LSR */
4a08d54d3a Mart*2128 bet = BV(i,j,bi,bj)-AV(i,j,bi,bj)*CVV(i,jM)
2129 CVV(i,j) = CVV(i,j)/bet
2130 VRT(i,j) = (VRT(i,j)-AV(i,j,bi,bj)*VRT(i,jM))/bet
5fb12930f9 Mart*2131 ENDDO
8d9be9cde5 Mart*2132 #ifdef SEAICE_VECTORIZE_LSR
2133 ENDDO
2134
032b010f0c Mart*2135
4a08d54d3a Mart*2136 DO j=jMin,jMax-1
2137 jM=sNy-j
2138 DO i=iMinLoc,iMax,iStep
032b010f0c Mart*2139 #else
2140
4a08d54d3a Mart*2141 DO j=jMin,jMax-1
2142 jM=sNy-j
032b010f0c Mart*2143 #endif /* SEAICE_VECTORIZE_LSR */
4a08d54d3a Mart*2144 VRT(i,jM)=VRT(i,jM)-CVV(i,jM)*VRT(i,jM+1)
5fb12930f9 Mart*2145 ENDDO
032b010f0c Mart*2146 #ifdef SEAICE_VECTORIZE_LSR
2147 ENDDO
2148
2149
4a08d54d3a Mart*2150 DO j=jMin,jMax
2151 DO i=iMinLoc,iMax,iStep
032b010f0c Mart*2152 #else
4a08d54d3a Mart*2153 DO j=jMin,jMax
032b010f0c Mart*2154 #endif /* SEAICE_VECTORIZE_LSR */
4a08d54d3a Mart*2155 vIce(i,j,bi,bj)=vTmp(i,j,bi,bj)
2156 & +WFAV*(VRT(i,j)-vTmp(i,j,bi,bj))
5fb12930f9 Mart*2157 ENDDO
2158 ENDDO
12dace3fd6 Mart*2159 #ifdef SEAICE_LSR_ZEBRA
8d9be9cde5 Mart*2160 ENDDO
12dace3fd6 Mart*2161 #endif /* SEAICE_LSR_ZEBRA */
5fb12930f9 Mart*2162
b7c2bb63b7 Mart*2163 #endif /* SEAICE_CGRID */
2164
2165 RETURN
2166 END