C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_lsr.F,v 1.124 2016/05/17 15:26:46 mlosch Exp $ C $Name: $ #include "SEAICE_OPTIONS.h" #ifdef ALLOW_OBCS # include "OBCS_OPTIONS.h" #else # define OBCS_UVICE_OLD #endif #ifdef ALLOW_AUTODIFF # include "AUTODIFF_OPTIONS.h" #endif C-- File seaice_lsr.F: seaice LSR dynamical solver S/R: C-- Contents C-- o SEAICE_LSR C-- o SEAICE_RESIDUAL C-- o SEAICE_LSR_CALC_COEFFS C-- o SEAICE_LSR_RHSU C-- o SEAICE_LSR_RHSV C-- o SEAICE_LSR_TRIDIAGU C-- o SEAICE_LSR_TRIDIAGV CBOP C !ROUTINE: SEAICE_LSR C !INTERFACE: SUBROUTINE SEAICE_LSR( myTime, myIter, myThid ) C !DESCRIPTION: \bv C *==========================================================* C | S/R SEAICE_LSR C | o Solve ice momentum equation with an LSR dynamics solver C | (see Zhang and Hibler, JGR, 102, 8691-8702, 1997 C | and Zhang and Rothrock, MWR, 131, 845- 861, 2003) C | Written by Jinlun Zhang, PSC/UW, Feb-2001 C | zhang@apl.washington.edu C *==========================================================* C | C-grid version by Martin Losch C | Since 2009/03/18: finite-Volume discretization of C | stress divergence that includes all metric terms C *==========================================================* C \ev C !USES: IMPLICIT NONE C === Global variables === #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "DYNVARS.h" #include "GRID.h" #include "SEAICE_SIZE.h" #include "SEAICE_PARAMS.h" #include "SEAICE.h" #ifdef ALLOW_AUTODIFF_TAMC # include "AUTODIFF_PARAMS.h" # include "tamc.h" #endif C !INPUT/OUTPUT PARAMETERS: C === Routine arguments === C myTime :: Simulation time C myIter :: Simulation timestep number C myThid :: my Thread Id. number _RL myTime INTEGER myIter INTEGER myThid CEOP #ifdef SEAICE_CGRID #ifdef SEAICE_ALLOW_DYNAMICS C !FUNCTIONS: LOGICAL DIFFERENT_MULTIPLE EXTERNAL DIFFERENT_MULTIPLE C !LOCAL VARIABLES: C === Local variables === C i,j,bi,bj :: Loop counters INTEGER ipass, ipass0 INTEGER i, j, m, bi, bj INTEGER iMin, iMax, jMin, jMax INTEGER ICOUNT1, ICOUNT2, linearIterLoc, nonLinIterLoc INTEGER kSrf INTEGER SOLV_MAX_TMP LOGICAL doIterate4u, doIterate4v CHARACTER*(MAX_LEN_MBUF) msgBuf _RL WFAU, WFAV, WFAU2, WFAV2 _RL S1, S2, S1A, S2A _RL eplus, eminus _RL recip_deltaT C diagonals of coefficient matrices _RL AU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL BU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL CU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL AV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL BV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL CV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C RHS of TriDiag Linear system (both directions => 5 coef A,B,C & Rt1,2) _RL rhsU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL rhsV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) #ifdef SEAICE_GLOBAL_3DIAG_SOLVER C RHS of TriDiag Linear system (1 direction only => 3 coef A,B,C) _RL rhsUx(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL rhsVy(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) INTEGER iSubIter INTEGER errCode #endif C coefficients for lateral points, u(j+/-1) _RL uRt1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL uRt2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C coefficients for lateral points, v(i+/-1) _RL vRt1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL vRt2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C abbreviations _RL etaPlusZeta (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL zetaMinusEta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C symmetric drag coefficient _RL dragSym (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C contribution of sigma on righ hand side _RL sig11(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL sig22(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL sig12(1-OLx:sNx+OLx,1-OLy:sNy+OLy) C intermediate time level centered (C) between n and n-1 _RL uIceC(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL vIceC(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C auxillary fields _RL uTmp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL vTmp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C fractional area at velocity points _RL areaW(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL areaS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) #ifdef ALLOW_AUTODIFF_TAMC INTEGER itmpkey, itmpkey2, itmpkey3 #endif _RL COSWAT _RS SINWAT _RL UERR #ifdef SEAICE_ALLOW_CHECK_LSR_CONVERGENCE _RL resnorm, EKnorm, counter #endif /* SEAICE_ALLOW_CHECK_LSR_CONVERGENCE */ LOGICAL printResidual _RL residUini, residVini, residUend, residVend #ifdef SEAICE_ALLOW_FREEDRIFT _RL uRes (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL vRes (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL errIni, errFD, errSum _RL residU_fd, residV_fd _RL residUmix, residVmix #endif /* SEAICE_ALLOW_FREEDRIFT */ C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| printResidual = debugLevel.GE.debLevA & .AND. DIFFERENT_MULTIPLE( SEAICE_monFreq, myTime, deltaTClock ) C extra overlap for (restricted) additive Schwarz method jMin = 1 - SEAICE_OLy jMax = sNy + SEAICE_OLy iMin = 1 - SEAICE_OLx iMax = sNx + SEAICE_OLx C linearIterLoc = SEAICElinearIterMax nonLinIterLoc = SEAICEnonLinIterMax IF ( SEAICEusePicardAsPrecon ) THEN linearIterLoc = SEAICEpreconlinIter nonLinIterLoc = SEAICEpreconNL_Iter ENDIF C recip_deltaT = 1. _d 0 / SEAICE_deltaTdyn #ifdef ALLOW_AUTODIFF_TAMC cph break artificial dependencies DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx deltaC (i,j,bi,bj) = 0. _d 0 press (i,j,bi,bj) = 0. _d 0 zeta (i,j,bi,bj) = 0. _d 0 zetaZ (i,j,bi,bj) = 0. _d 0 eta (i,j,bi,bj) = 0. _d 0 etaZ (i,j,bi,bj) = 0. _d 0 uIceC (i,j,bi,bj) = 0. _d 0 vIceC (i,j,bi,bj) = 0. _d 0 uIceNm1(i,j,bi,bj) = 0. _d 0 vIceNm1(i,j,bi,bj) = 0. _d 0 ENDDO ENDDO ENDDO ENDDO #endif C initialise fractional areas at velocity points DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx areaW(I,J,bi,bj) = 1. _d 0 areaS(I,J,bi,bj) = 1. _d 0 ENDDO ENDDO IF ( SEAICEscaleSurfStress ) THEN DO J=jMin,jMax DO I=iMin,iMax areaW(I,J,bi,bj) = & 0.5 _d 0*(AREA(I,J,bi,bj)+AREA(I-1,J,bi,bj)) areaS(I,J,bi,bj) = & 0.5 _d 0*(AREA(I,J,bi,bj)+AREA(I,J-1,bi,bj)) ENDDO ENDDO ENDIF ENDDO ENDDO ipass0 = 2 C modify time stepping a little if more than the default C modified Eulerian time step is used. IF ( nonLinIterLoc.GT.2 ) ipass0 = 1 #ifdef ALLOW_AUTODIFF_TAMC DO ipass=1,MPSEUDOTIMESTEPS itmpkey = (ikey_dynamics-1)*MPSEUDOTIMESTEPS + ipass #else DO ipass=1,nonLinIterLoc #endif #ifdef ALLOW_AUTODIFF_TAMC # ifdef SEAICE_ALLOW_FREEDRIFT CADJ STORE uice_fd, vice_fd = comlev1_dynsol, kind=isbyte, CADJ & key = itmpkey CADJ STORE ures,vres = comlev1_dynsol, kind=isbyte, CADJ & key = itmpkey # endif CADJ STORE utmp,vtmp = comlev1_dynsol, kind=isbyte, CADJ & key = itmpkey #endif /* ALLOW_AUTODIFF_TAMC */ IF ( ipass .LE. nonLinIterLoc ) THEN #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE uice = comlev1_dynsol, kind=isbyte, CADJ & key = itmpkey CADJ STORE vice = comlev1_dynsol, kind=isbyte, CADJ & key = itmpkey CADJ STORE uicenm1 = comlev1_dynsol, kind=isbyte, CADJ & key = itmpkey CADJ STORE vicenm1 = comlev1_dynsol, kind=isbyte, CADJ & key = itmpkey #endif /* ALLOW_AUTODIFF_TAMC */ C surface level kSrf = 1 C-- introduce turning angles SINWAT=SIN(SEAICE_waterTurnAngle*deg2rad) COSWAT=COS(SEAICE_waterTurnAngle*deg2rad) C SET SOME VALUES WFAU=SEAICE_LSRrelaxU WFAV=SEAICE_LSRrelaxV WFAU2=ZERO WFAV2=ZERO S1 = 0. S2 = 0. S1A=0.80 _d 0 S2A=0.80 _d 0 IF ( ipass .EQ. 1 ) THEN C NOW DO PREDICTOR TIME STEP DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx uIceNm1(I,J,bi,bj) = uIce(I,J,bi,bj) vIceNm1(I,J,bi,bj) = vIce(I,J,bi,bj) uIceC (I,J,bi,bj) = uIce(I,J,bi,bj) vIceC (I,J,bi,bj) = vIce(I,J,bi,bj) ENDDO ENDDO ENDDO ENDDO ELSE C NOW DO MODIFIED EULER STEP DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx uIce (I,J,bi,bj)=HALF*(uIce(I,J,bi,bj)+uIceNm1(i,j,bi,bj)) vIce (I,J,bi,bj)=HALF*(vIce(I,J,bi,bj)+vIceNm1(i,j,bi,bj)) uIceC(I,J,bi,bj)=uIce(I,J,bi,bj) vIceC(I,J,bi,bj)=vIce(I,J,bi,bj) ENDDO ENDDO ENDDO ENDDO ENDIF IF ( ipass .GT. ipass0 ) THEN C for additional (pseudo-time)steps update u/vIceNm1 DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx uIceNm1(i,j,bi,bj)=uIce(I,J,bi,bj) vIceNm1(i,j,bi,bj)=vIce(I,J,bi,bj) ENDDO ENDDO ENDDO ENDDO ENDIF #ifdef ALLOW_AUTODIFF_TAMC C That is an important one! Note, that C * lsr is called multiple times (twice by default), thus the C ipass index is part of itmpkey C * this storing is still outside the lsr iteration loop CADJ STORE uice = comlev1_dynsol, kind=isbyte, key = itmpkey CADJ STORE vice = comlev1_dynsol, kind=isbyte, key = itmpkey #endif /* ALLOW_AUTODIFF_TAMC */ CALL SEAICE_CALC_STRAINRATES( I uIceC, vIceC, O e11, e22, e12, I ipass, myTime, myIter, myThid ) CALL SEAICE_CALC_VISCOSITIES( I e11, e22, e12, zMin, zMax, hEffM, press0, tensileStrFac, O eta, etaZ, zeta, zetaZ, press, deltaC, I ipass, myTime, myIter, myThid ) CALL SEAICE_OCEANDRAG_COEFFS( I uIceC, vIceC, O DWATN, I ipass, myTime, myIter, myThid ) #ifdef SEAICE_ALLOW_BOTTOMDRAG IF (SEAICEbasalDragK2.GT.0. _d 0) CALL SEAICE_BOTTOMDRAG_COEFFS( I uIceC, vIceC, #ifdef SEAICE_ITD I HEFFITD, AREAITD, AREA, #else I HEFF, AREA, #endif O CbotC, I ipass, myTime, myIter, myThid ) #endif /* SEAICE_ALLOW_BOTTOMDRAG */ C C some abbreviations C DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO J=jMin-1,jMax DO I=iMin-1,iMax etaPlusZeta (I,J,bi,bj) = ETA (I,J,bi,bj)+ZETA(I,J,bi,bj) zetaMinusEta(I,J,bi,bj) = ZETA(I,J,bi,bj)-ETA (I,J,bi,bj) ENDDO ENDDO DO J=1-OLy,sNy+OLy DO I=1-OLx,sNx+OLx dragSym(I,J,bi,bj) = DWATN(I,J,bi,bj)*COSWAT #ifdef SEAICE_ALLOW_BOTTOMDRAG & +CbotC(I,J,bi,bj) #endif /* SEAICE_ALLOW_BOTTOMDRAG */ ENDDO ENDDO ENDDO ENDDO C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| C C Set up of right hand sides of momentum equations starts here C DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO J=jMin,jMax DO I=iMin,iMax C set up anti symmetric drag force and add in current force C ( remember to average to correct velocity points ) FORCEX(I,J,bi,bj)=FORCEX0(I,J,bi,bj)+ & ( 0.5 _d 0 * ( DWATN(I,J,bi,bj)+DWATN(I-1,J,bi,bj) ) * & COSWAT * uVel(I,J,kSrf,bi,bj) & - SIGN(SINWAT, _fCori(I,J,bi,bj))* 0.5 _d 0 * & ( DWATN(I ,J,bi,bj) * 0.5 _d 0 * & (vVel(I ,J ,kSrf,bi,bj)-vIceC(I ,J ,bi,bj) & +vVel(I ,J+1,kSrf,bi,bj)-vIceC(I ,J+1,bi,bj)) & + DWATN(I-1,J,bi,bj) * 0.5 _d 0 * & (vVel(I-1,J ,kSrf,bi,bj)-vIceC(I-1,J ,bi,bj) & +vVel(I-1,J+1,kSrf,bi,bj)-vIceC(I-1,J+1,bi,bj)) & ) ) * areaW(I,J,bi,bj) FORCEY(I,J,bi,bj)=FORCEY0(I,J,bi,bj)+ & ( 0.5 _d 0 * ( DWATN(I,J,bi,bj)+DWATN(I,J-1,bi,bj) ) * & COSWAT * vVel(I,J,kSrf,bi,bj) & + SIGN(SINWAT, _fCori(I,J,bi,bj)) * 0.5 _d 0 * & ( DWATN(I,J ,bi,bj) * 0.5 _d 0 * & (uVel(I ,J ,kSrf,bi,bj)-uIceC(I ,J ,bi,bj) & +uVel(I+1,J ,kSrf,bi,bj)-uIceC(I+1,J ,bi,bj)) & + DWATN(I,J-1,bi,bj) * 0.5 _d 0 * & (uVel(I ,J-1,kSrf,bi,bj)-uIceC(I ,J-1,bi,bj) & +uVel(I+1,J-1,kSrf,bi,bj)-uIceC(I+1,J-1,bi,bj)) & ) ) * areaS(I,J,bi,bj) ENDDO ENDDO DO j=jMin,jMax DO i=iMin,iMax C- add Coriolis term FORCEX(I,J,bi,bj) = FORCEX(I,J,bi,bj) + HALF* & ( seaiceMassC(I ,J,bi,bj) * _fCori(I ,J,bi,bj) & *0.5 _d 0*( vIceC( i ,j,bi,bj)+vIceC( i ,j+1,bi,bj) ) & + seaiceMassC(I-1,J,bi,bj) * _fCori(I-1,J,bi,bj) & *0.5 _d 0*( vIceC(i-1,j,bi,bj)+vIceC(i-1,j+1,bi,bj) ) ) FORCEY(I,J,bi,bj) = FORCEY(I,J,bi,bj) - HALF* & ( seaiceMassC(I,J ,bi,bj) * _fCori(I,J ,bi,bj) & *0.5 _d 0*( uIceC(i ,j ,bi,bj)+uIceC(i+1, j,bi,bj) ) & + seaiceMassC(I,J-1,bi,bj) * _fCori(I,J-1,bi,bj) & *0.5 _d 0*( uIceC(i ,j-1,bi,bj)+uIceC(i+1,j-1,bi,bj) ) ) ENDDO ENDDO C this is the rhs contribution of the time derivative DO j=jMin,jMax DO i=iMin,iMax FORCEX(I,J,bi,bj)=FORCEX(I,J,bi,bj) & +seaiceMassU(I,J,bi,bj)*recip_deltaT & *uIceNm1(i,j,bi,bj) FORCEY(I,J,bi,bj)=FORCEY(I,J,bi,bj) & +seaiceMassV(I,J,bi,bj)*recip_deltaT & *vIceNm1(i,j,bi,bj) FORCEX(I,J,bi,bj)=FORCEX(I,J,bi,bj)* seaiceMaskU(I,J,bi,bj) FORCEY(I,J,bi,bj)=FORCEY(I,J,bi,bj)* seaiceMaskV(I,J,bi,bj) ENDDO ENDDO C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| C C u-equation C DO j=jMin,jMax DO i=iMin,iMax rhsU(I,J,bi,bj) = FORCEX(I,J,bi,bj) ENDDO ENDDO CALL SEAICE_LSR_RHSU( I zetaMinusEta, etaPlusZeta, etaZ, zetaZ, press, I uIceC, vIceC, U rhsU, I iMin, iMax, jMin, jMax, bi, bj, myThid ) C C v-equation C DO j=jMin,jMax DO i=iMin,iMax rhsV(I,J,bi,bj) = FORCEY(I,J,bi,bj) ENDDO ENDDO CALL SEAICE_LSR_RHSV( I zetaMinusEta, etaPlusZeta, etaZ, zetaZ, press, I uIceC, vIceC, U rhsV, I iMin, iMax, jMin, jMax, bi, bj, myThid ) ENDDO ENDDO C C Set up of right hand sides of momentum equations ends here C C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| C C calculate coefficients of tridiagonal matrices for both u- and C v-equations C CALL SEAICE_LSR_CALC_COEFFS( I etaPlusZeta, zetaMinusEta, etaZ, zetaZ, dragSym, O AU, BU, CU, AV, BV, CV, uRt1, uRt2, vRt1, vRt2, I iMin, iMax, jMin, jMax, myTime, myIter, myThid ) #ifndef OBCS_UVICE_OLD C-- prevent tri-diagonal solver from modifying OB values: DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO J=jMin,jMax DO I=iMin,iMax IF ( maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj) .EQ. 0. ) THEN AU(I,J,bi,bj) = ZERO BU(I,J,bi,bj) = ONE CU(I,J,bi,bj) = ZERO uRt1(I,J,bi,bj) = ZERO uRt2(I,J,bi,bj) = ZERO rhsU(I,J,bi,bj) = uIce(I,J,bi,bj) ENDIF IF ( maskInC(i,j,bi,bj)*maskInC(i,j-1,bi,bj) .EQ. 0. ) THEN AV(I,J,bi,bj) = ZERO BV(I,J,bi,bj) = ONE CV(I,J,bi,bj) = ZERO vRt1(I,J,bi,bj) = ZERO vRt2(I,J,bi,bj) = ZERO rhsV(I,J,bi,bj) = vIce(I,J,bi,bj) ENDIF ENDDO ENDDO ENDDO ENDDO #endif /* OBCS_UVICE_OLD */ C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| #ifdef ALLOW_DEBUG IF ( debugLevel .GE. debLevD ) THEN WRITE(msgBuf,'(A,I3,A)') & 'Uice pre iter (SEAICE_LSR', MOD(ipass,1000), ')' CALL DEBUG_STATS_RL( 1, UICE, msgBuf, myThid ) WRITE(msgBuf,'(A,I3,A)') & 'Vice pre iter (SEAICE_LSR', MOD(ipass,1000), ')' CALL DEBUG_STATS_RL( 1, VICE, msgBuf, myThid ) ENDIF #endif /* ALLOW_DEBUG */ C-- Calculate initial residual of the linearised system IF ( printResidual .OR. LSR_mixIniGuess.GE.1 ) & CALL SEAICE_RESIDUAL( I rhsU, rhsV, uRt1, uRt2, vRt1, vRt2, I AU, BU, CU, AV, BV, CV, uIce, vIce, O residUini, residVini, uTmp, vTmp, I printResidual, myIter, myThid ) #ifdef SEAICE_ALLOW_FREEDRIFT IF ( printResidual .OR. LSR_mixIniGuess.GE.1 ) THEN IF ( LSR_mixIniGuess.GE.1 ) THEN DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=jMin,jMax DO i=iMin,iMax uIce_fd(i,j,bi,bj) = FORCEX(i,j,bi,bj) & / ( 1. _d 0 - seaiceMaskU(i,j,bi,bj) & + seaiceMassU(i,j,bi,bj)*recip_deltaT & + HALF*( dragSym(i,j,bi,bj) + dragSym(i-1,j,bi,bj) ) & * areaW(i,j,bi,bj) & ) vIce_fd(i,j,bi,bj) = FORCEY(i,j,bi,bj) & / ( 1. _d 0 - seaiceMaskV(i,j,bi,bj) & + seaiceMassV(i,j,bi,bj)*recip_deltaT & + HALF*( dragSym(i,j,bi,bj) + dragSym(i,j-1,bi,bj) ) & * areaS(i,j,bi,bj) & ) ENDDO ENDDO ENDDO ENDDO CALL EXCH_UV_XY_RL( uIce_fd, vIce_fd, .TRUE., myThid ) ENDIF IF ( LSR_mixIniGuess.GE.0 ) & CALL SEAICE_RESIDUAL( I rhsU, rhsV, uRt1, uRt2, vRt1, vRt2, I AU, BU, CU, AV, BV, CV, uIce_fd, vIce_fd, O residU_fd, residV_fd, uRes, vRes, I printResidual, myIter, myThid ) ENDIF IF ( LSR_mixIniGuess.GE.2 ) THEN # ifdef ALLOW_AUTODIFF_TAMC CADJ STORE uice, vice = comlev1_dynsol, kind=isbyte, CADJ & key = itmpkey CADJ STORE utmp, vtmp = comlev1_dynsol, kind=isbyte, CADJ & key = itmpkey # endif /* ALLOW_AUTODIFF_TAMC */ DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=jMin,jMax DO i=iMin,iMax errIni = uTmp(i,j,bi,bj)*uTmp(i,j,bi,bj) errFD = uRes(i,j,bi,bj)*uRes(i,j,bi,bj) IF ( LSR_mixIniGuess.GE.4 ) THEN errIni = errIni*errIni errFD = errFD *errFD ENDIF errSum = ( errIni + errFD ) & *maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj) IF ( errSum.GT.0. ) THEN uIce(i,j,bi,bj) = ( errFD *uIce(i,j,bi,bj) & + errIni*uIce_fd(i,j,bi,bj) & )/errSum ENDIF errIni = vTmp(i,j,bi,bj)*vTmp(i,j,bi,bj) errFD = vRes(i,j,bi,bj)*vRes(i,j,bi,bj) IF ( LSR_mixIniGuess.GE.4 ) THEN errIni = errIni*errIni errFD = errFD *errFD ENDIF errSum = ( errIni + errFD ) & *maskInC(i,j,bi,bj)*maskInC(i,j-1,bi,bj) IF ( errSum.GT.0. ) THEN vIce(i,j,bi,bj) = ( errFD *vIce(i,j,bi,bj) & + errIni*vIce_fd(i,j,bi,bj) & )/errSum ENDIF ENDDO ENDDO ENDDO ENDDO CALL EXCH_UV_XY_RL( uIce, vIce, .TRUE., myThid ) #ifndef ALLOW_AUTODIFF_TAMC C Here residU/Vmix and u/vTmp are computed but never used for anything C but reporting to STDOUT. Still TAF thinks this requires recomputations C so we hide this part of the code from TAF. IF ( printResidual ) THEN CALL SEAICE_RESIDUAL( I rhsU, rhsV, uRt1, uRt2, vRt1, vRt2, I AU, BU, CU, AV, BV, CV, uIce, vIce, O residUmix, residVmix, uTmp, vTmp, I printResidual, myIter, myThid ) ENDIF #endif /* ALLOW_AUTODIFF_TAMC */ ENDIF #endif /* SEAICE_ALLOW_FREEDRIFT */ C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| C NOW DO ITERATION cph--- iteration starts here cph--- need to kick out goto doIterate4u = .TRUE. doIterate4v = .TRUE. C ITERATION START ----------------------------------------------------- C Ideally, we would like to use the loop-directive for this C iterative loop, but (unsurprisingly) it does not seem to work for C this convoluted case of mixed iterations, so it is commented out CML#ifdef ALLOW_AUTODIFF_TAMC CMLCADJ LOOP = iteration uice, vice = comlev1_lsr CML#endif /* ALLOW_AUTODIFF_TAMC */ #if (defined (ALLOW_AUTODIFF_TAMC) && defined (SEAICE_LSR_ADJOINT_ITER)) IF ( inAdMode ) THEN SOLV_MAX_TMP = SOLV_MAX_FIXED ELSE SOLV_MAX_TMP = linearIterLoc ENDIF #else SOLV_MAX_TMP = linearIterLoc #endif ICOUNT1 = SOLV_MAX_TMP ICOUNT2 = SOLV_MAX_TMP DO m = 1, SOLV_MAX_TMP #ifdef ALLOW_AUTODIFF_TAMC # ifdef SEAICE_LSR_ADJOINT_ITER itmpkey2 = (itmpkey-1)*SOLV_MAX_FIXED + MIN(m,SOLV_MAX_FIXED) # else itmpkey2 = itmpkey # endif /* SEAICE_LSR_ADJOINT_ITER */ C Storing these scalars and logicals is necessary to do an exact C backward iteration CADJ STORE wfau, wfav = comlev1_lsr, kind=isbyte, key = itmpkey2 CADJ STORE doIterate4u, doIterate4v = comlev1_lsr, key = itmpkey2 #endif /* ALLOW_AUTODIFF_TAMC */ IF ( doIterate4u .OR. doIterate4v ) THEN IF ( useCubedSphereExchange ) THEN doIterate4u = .TRUE. doIterate4v = .TRUE. ENDIF # ifdef ALLOW_AUTODIFF_TAMC C Note that itmpkey2 can get very large when SEAICE_LSR_ADJOINT_ITER C is defined for an accurate adjoint, otherwise it is the same for C each m, so that the tapes get overwritten for each iterate and C the adjoint is necessarily wrong. CADJ STORE uice, vice = comlev1_lsr, kind=isbyte, key = itmpkey2 # endif /* ALLOW_AUTODIFF_TAMC */ #ifdef SEAICE_GLOBAL_3DIAG_SOLVER IF ( SEAICEuseMultiTileSolver ) THEN DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) IF ( doIterate4u ) THEN DO j=jMin,jMax DO i=iMin,iMax uTmp(i,j,bi,bj) = uIce(i,j,bi,bj) rhsUx(i,j,bi,bj) = ( rhsU(i,j,bi,bj) & + uRt1(i,j,bi,bj)*uIce(i,j-1,bi,bj) & + uRt2(i,j,bi,bj)*uIce(i,j+1,bi,bj) & )*seaiceMaskU(i,j,bi,bj) ENDDO ENDDO ENDIF IF ( doIterate4v ) THEN DO j=jMin,jMax DO i=iMin,iMax vTmp(i,j,bi,bj) = vIce(i,j,bi,bj) rhsVy(i,j,bi,bj) = ( rhsV(i,j,bi,bj) & + vRt1(i,j,bi,bj)*vIce(i-1,J,bi,bj) & + vRt2(i,j,bi,bj)*vIce(i+1,J,bi,bj) & )*seaiceMaskU(i,j,bi,bj) ENDDO ENDDO ENDIF C end bi,bj-loops ENDDO ENDDO iSubIter = SOLV_MAX_TMP*(ipass-1) + m CALL SOLVE_UV_TRIDIAGO( I 1, OLx, doIterate4u, doIterate4v, I AU, BU, CU, rhsUx, I AV, BV, CV, rhsVy, O uIce, vIce, errCode, I iSubIter, myIter, myThid ) DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) IF ( doIterate4u ) THEN DO j=jMin,jMax DO i=iMin,iMax uIce(i,j,bi,bj) = uTmp(i,j,bi,bj) & + WFAU*( uIce(i,j,bi,bj)-uTmp(i,j,bi,bj) ) ENDDO ENDDO ENDIF IF ( doIterate4v ) THEN DO j=jMin,jMax DO i=iMin,iMax vIce(i,j,bi,bj) = vTmp(i,j,bi,bj) & + WFAV*( vIce(i,j,bi,bj)-vTmp(i,j,bi,bj) ) ENDDO ENDDO ENDIF C end bi,bj-loops ENDDO ENDDO ELSE #endif /* SEAICE_GLOBAL_3DIAG_SOLVER */ DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) #ifdef ALLOW_AUTODIFF_TAMC act1 = bi - myBxLo(myThid) max1 = myBxHi(myThid) - myBxLo(myThid) + 1 act2 = bj - myByLo(myThid) max2 = myByHi(myThid) - myByLo(myThid) + 1 act3 = myThid - 1 max3 = nTx*nTy act4 = itmpkey2 - 1 itmpkey3 = (act1 + 1) + act2*max1 & + act3*max1*max2 & + act4*max1*max2*max3 CADJ STORE uice(:,:,bi,bj) = comlev1_bibj_lsr,kind=isbyte,key=itmpkey3 CADJ STORE vice(:,:,bi,bj) = comlev1_bibj_lsr,kind=isbyte,key=itmpkey3 #endif /* ALLOW_AUTODIFF_TAMC */ C-jmc: get less TAF warnings when always (no if doIterate) saving uIce,vIce: C save u/vIce prior to iteration DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx uTmp(I,J,bi,bj)=uIce(I,J,bi,bj) vTmp(I,J,bi,bj)=vIce(I,J,bi,bj) ENDDO ENDDO IF ( doIterate4u ) THEN C Solve for uIce : CALL SEAICE_LSR_TRIDIAGU( I AU, BU, CU, uRt1, uRt2, rhsU, uTmp, seaiceMaskU, WFAU, U uIce, I iMin, iMax, jMin, jMax, bi, bj, myTime, myIter, myThid ) ENDIF IF ( doIterate4v ) THEN C Solve for vIce CALL SEAICE_LSR_TRIDIAGV( I AV, BV, CV, vRt1, vRt2, rhsV, vTmp, seaiceMaskV, WFAV, U vIce, I iMin, iMax, jMin, jMax, bi, bj, myTime, myIter, myThid ) ENDIF C end bi,bj-loops ENDDO ENDDO #ifdef SEAICE_GLOBAL_3DIAG_SOLVER ENDIF #endif /* SEAICE_GLOBAL_3DIAG_SOLVER */ IF ( doIterate4u.AND.MOD(m,SOLV_NCHECK).EQ.0) THEN S1=ZERO DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO J=1,sNy DO I=1,sNx UERR=(uIce(I,J,bi,bj)-uTmp(I,J,bi,bj)) & * seaiceMaskU(I,J,bi,bj) S1=MAX(ABS(UERR),S1) ENDDO ENDDO ENDDO ENDDO _GLOBAL_MAX_RL( S1, myThid ) c WRITE(standardMessageUnit,'(A,2I6,1P4E16.9)') c & ' U iters,error,WF = ',ipass,M,S1,S1A,WFAU C SAFEGUARD AGAINST BAD FORCING ETC IF(m.GT.1.AND.S1.GT.S1A) WFAU=WFAU2 S1A=S1 IF(S1.LT.LSR_ERROR) THEN ICOUNT1=m doIterate4u = .FALSE. ENDIF ENDIF IF ( doIterate4v.AND.MOD(m,SOLV_NCHECK).EQ.0) THEN S2=ZERO DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO J=1,sNy DO I=1,sNx UERR=(vIce(I,J,bi,bj)-vTmp(I,J,bi,bj)) & * seaiceMaskV(I,J,bi,bj) S2=MAX(ABS(UERR),S2) ENDDO ENDDO ENDDO ENDDO _GLOBAL_MAX_RL( S2, myThid ) C SAFEGUARD AGAINST BAD FORCING ETC IF(m.GT.1.AND.S2.GT.S2A) WFAV=WFAV2 S2A=S2 IF(S2.LT.LSR_ERROR) THEN ICOUNT2=m doIterate4v = .FALSE. ENDIF ENDIF CALL EXCH_UV_XY_RL( uIce, vIce,.TRUE.,myThid) C-- end doIterate4u or doIterate4v ENDIF ENDDO C ITERATION END ----------------------------------------------------- C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| #ifndef ALLOW_AUTODIFF_TAMC C Here residU/Vend and u/vTmp are computed but never used for anything C but reporting to STDOUT. Still TAF thinks this requires recomputations C so we hide this part of the code from TAF. IF ( printResidual ) THEN C-- Calculate final residual of the linearised system CALL SEAICE_RESIDUAL( I rhsU, rhsV, uRt1, uRt2, vRt1, vRt2, I AU, BU, CU, AV, BV, CV, uIce, vIce, O residUend, residVend, uTmp, vTmp, I printResidual, myIter, myThid ) _BEGIN_MASTER( myThid ) WRITE(standardMessageUnit,'(A,1X,I4,1P3E16.8)') & ' SEAICE_LSR: Residual Initial ipass,Uice,Vice=', & ipass, residUini, residVini #ifdef SEAICE_ALLOW_FREEDRIFT IF ( LSR_mixIniGuess.GE.0 ) & WRITE(standardMessageUnit,'(A,1P3E16.8)') & ' SEAICE_LSR: Residual FrDrift U_fd,V_fd=', & residU_fd, residV_fd IF ( LSR_mixIniGuess.GE.2 ) & WRITE(standardMessageUnit,'(A,1P3E16.8)') & ' SEAICE_LSR: Residual Mix-ini Uice,Vice=', & residUmix, residVmix #endif /* SEAICE_ALLOW_FREEDRIFT */ WRITE(standardMessageUnit,'(A,I4,A,I6,1P2E16.8)') & ' SEAICE_LSR (ipass=',ipass,') iters,dU,Resid=', & ICOUNT1, S1, residUend WRITE(standardMessageUnit,'(A,I4,A,I6,1P2E16.8)') & ' SEAICE_LSR (ipass=',ipass,') iters,dV,Resid=', & ICOUNT2, S2, residVend _END_MASTER( myThid ) ENDIF #ifdef ALLOW_DEBUG IF ( debugLevel .GE. debLevD ) THEN WRITE(msgBuf,'(A,I3,A)') & 'Uice post iter (SEAICE_LSR', MOD(ipass,1000), ')' CALL DEBUG_STATS_RL( 1, UICE, msgBuf, myThid ) WRITE(msgBuf,'(A,I3,A)') & 'Vice post iter (SEAICE_LSR', MOD(ipass,1000), ')' CALL DEBUG_STATS_RL( 1, VICE, msgBuf, myThid ) ENDIF #endif /* ALLOW_DEBUG */ IF ( doIterate4u .OR. doIterate4v ) THEN WRITE(msgBuf,'(2A,I10,A,I4,A)') '** WARNING ** SEAICE_LSR ', & '(it=', myIter, ',', ipass,') did not converge :' CALL PRINT_MESSAGE( msgBuf, errorMessageUnit, & SQUEEZE_RIGHT, myThid ) WRITE(msgBuf,'(2(A,I6,0PF6.3,1PE14.6))') & ' nIt,wFU,dU=', ICOUNT1, WFAU, S1, & ' ; nIt,wFV,dV=', ICOUNT2, WFAV, S2 CALL PRINT_MESSAGE( msgBuf, errorMessageUnit, & SQUEEZE_RIGHT, myThid ) ENDIF #endif /* ALLOW_AUTODIFF_TAMC */ C APPLY MASKS DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO J=1-OLy,sNy+OLy DO I=1-OLx,sNx+OLx uIce(I,J,bi,bj)=uIce(I,J,bi,bj)* seaiceMaskU(I,J,bi,bj) vIce(I,J,bi,bj)=vIce(I,J,bi,bj)* seaiceMaskV(I,J,bi,bj) ENDDO ENDDO ENDDO ENDDO CML CALL EXCH_UV_XY_RL( uIce, vIce,.TRUE.,myThid) #ifdef SEAICE_ALLOW_CHECK_LSR_CONVERGENCE IF ( debugLevel .GE. debLevC ) THEN resnorm = 0. _d 0 EKnorm = 0. _d 0 counter = 0. _d 0 DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) C-- Compute residual norm and print DO j=1,sNy DO i=1,sNx resnorm = resnorm + 0.5 _d 0 * & ( ( (uIceNm1(i,j,bi,bj)+uIceNm1(i+1,j,bi,bj)) & - (uIce(i,j,bi,bj)+uIce(i+1,j,bi,bj)) )**2 & + ( (vIceNm1(i,j,bi,bj)+vIceNm1(i,j+1,bi,bj)) & - (vIce(i,j,bi,bj)+vIce(i,j+1,bi,bj)) )**2 ) IF ( area(i,j,bi,bj) .gt. 0.5 _d 0 ) THEN EKnorm = EKnorm + 0.5 _d 0 * heff(i,j,bi,bj) * & ( ( (uIce(i,j,bi,bj)+uIce(i+1,j,bi,bj)) )**2 & + ( (vIce(i,j,bi,bj)+vIce(i,j+1,bi,bj)) )**2 ) counter = counter + 1. _d 0 ENDIF ENDDO ENDDO ENDDO ENDDO _GLOBAL_SUM_RL( resnorm, myThid ) _GLOBAL_SUM_RL( EKnorm, myThid ) _GLOBAL_SUM_RL( counter, myThid ) IF ( counter .gt. 0. _d 0 ) EKnorm = EKnorm/counter _BEGIN_MASTER( myThid ) WRITE(standardMessageUnit,'(A,I7,1X,2E22.14)') & 'S/R SEAICE_LSR: IPSEUDO, RESNORM, EKNORM = ', & ipass, sqrt(resnorm), EKnorm _END_MASTER( myThid ) ENDIF #endif /* SEAICE_ALLOW_CHECK_LSR_CONVERGENCE */ ENDIF C end outer ipass pseudo-time stepping loop ENDDO IF ( useHB87StressCoupling ) THEN # ifdef ALLOW_AUTODIFF_TAMC C These directives do not remove any recomputation warnings but avoid C recomputing the entire LSR loop before evaluating this code block CADJ STORE uice, vice, eta, etaz, zeta = comlev1_dynsol, CADJ & kind=isbyte, key = ikey_dynamics # endif /* ALLOW_AUTODIFF_TAMC */ C compute the divergence of stress here to be used later C C compute strain rate from latest velocities CALL SEAICE_CALC_STRAINRATES( I uIce, vIce, O e11, e22, e12, I 3, myTime, myIter, myThid ) C compute internal stresses with updated ice velocities DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx sig11(I,J) = 0. _d 0 sig22(I,J) = 0. _d 0 sig12(I,J) = 0. _d 0 ENDDO ENDDO DO j=0,sNy DO i=0,sNx eplus = e11(I,J,bi,bj) + e22(I,J,bi,bj) eminus= e11(I,J,bi,bj) - e22(I,J,bi,bj) sig11(I,J) = zeta(I,J,bi,bj)*eplus + eta(I,J,bi,bj)*eminus & - 0.5 _d 0 * PRESS(I,J,bi,bj) sig22(I,J) = zeta(I,J,bi,bj)*eplus - eta(I,J,bi,bj)*eminus & - 0.5 _d 0 * PRESS(I,J,bi,bj) ENDDO ENDDO DO j=1,sNy+1 DO i=1,sNx+1 sig12(I,J) = 2. _d 0 * e12(I,J,bi,bj) * etaZ(I,J,bi,bj) ENDDO ENDDO C evaluate divergence of stress and apply to forcing DO J=1,sNy DO I=1,sNx stressDivergenceX(I,J,bi,bj) = & ( sig11(I ,J ) * _dyF(I ,J ,bi,bj) & - sig11(I-1,J ) * _dyF(I-1,J ,bi,bj) & + sig12(I ,J+1) * _dxV(I ,J+1,bi,bj) & - sig12(I ,J ) * _dxV(I ,J ,bi,bj) & ) * recip_rAw(I,J,bi,bj) stressDivergenceY(I,J,bi,bj) = & ( sig22(I ,J ) * _dxF(I ,J ,bi,bj) & - sig22(I ,J-1) * _dxF(I ,J-1,bi,bj) & + sig12(I+1,J ) * _dyU(I+1,J ,bi,bj) & - sig12(I ,J ) * _dyU(I ,J ,bi,bj) & ) * recip_rAs(I,J,bi,bj) ENDDO ENDDO ENDDO ENDDO C endif useHB87StressCoupling ENDIF #endif /* SEAICE_ALLOW_DYNAMICS */ #endif /* SEAICE_CGRID */ RETURN END C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP C !ROUTINE: SEAICE_RESIDUAL C !INTERFACE: SUBROUTINE SEAICE_RESIDUAL( I rhsU, rhsV, uRt1, uRt2, vRt1, vRt2, I AU, BU, CU, AV, BV, CV, uFld, vFld, O residU, residV, uRes, vRes, I calcMeanResid, myIter, myThid ) C !DESCRIPTION: \bv C *==========================================================* C | S/R SEAICE_RESIDUAL C | o Compute Linear solver residual C *==========================================================* C \ev C !USES: IMPLICIT NONE C === Global variables === #include "SIZE.h" #include "EEPARAMS.h" #include "GRID.h" #include "SEAICE_SIZE.h" #include "SEAICE.h" C !INPUT/OUTPUT PARAMETERS: C === Routine arguments === C calcMeanResid :: also compute mean residual (=RMS) C myIter :: Simulation timestep number C myThid :: my Thread Id. number C-- RHS _RL rhsU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL rhsV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C- coefficients for lateral points, u(j+/-1) _RL uRt1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL uRt2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C- coefficients for lateral points, v(i+/-1) _RL vRt1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL vRt2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C- diagonals of coefficient matrices _RL AU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL BU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL CU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL AV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL BV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL CV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C-- seaice velocity _RL uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C- residual (output) _RL residU, residV _RL uRes (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL vRes (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) LOGICAL calcMeanResid INTEGER myIter INTEGER myThid CEOP #ifdef SEAICE_CGRID #ifdef SEAICE_ALLOW_DYNAMICS C !LOCAL VARIABLES: C === Local variables === C i,j,bi,bj :: Loop counters INTEGER i, j, bi, bj C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| C-- Calculate residual of the linearised system residU = 0. residV = 0. DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1,sNy DO i=1,sNx uRes(i,j,bi,bj) = rhsU(i,j,bi,bj) & + uRt1(i,j,bi,bj)*uFld(i,j-1,bi,bj) & + uRt2(i,j,bi,bj)*uFld(i,j+1,bi,bj) & - ( AU(i,j,bi,bj)*uFld(i-1,j,bi,bj) & + BU(i,j,bi,bj)*uFld( i ,j,bi,bj) & + CU(i,j,bi,bj)*uFld(i+1,j,bi,bj) & ) vRes(i,j,bi,bj) = rhsV(i,j,bi,bj) & + vRt1(i,j,bi,bj)*vFld(i-1,j,bi,bj) & + vRt2(i,j,bi,bj)*vFld(i+1,j,bi,bj) & - ( AV(i,j,bi,bj)*vFld(i,j-1,bi,bj) & + BV(i,j,bi,bj)*vFld(i, j ,bi,bj) & + CV(i,j,bi,bj)*vFld(i,j+1,bi,bj) & ) ENDDO ENDDO IF ( calcMeanResid ) THEN DO j=1,sNy DO i=1,sNx residU = residU + uRes(i,j,bi,bj)*uRes(i,j,bi,bj) & *rAw(i,j,bi,bj)*maskInW(i,j,bi,bj) #ifndef OBCS_UVICE_OLD & *maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj) #endif residV = residV + vRes(i,j,bi,bj)*vRes(i,j,bi,bj) & *rAs(i,j,bi,bj)*maskInS(i,j,bi,bj) #ifndef OBCS_UVICE_OLD & *maskInC(i,j,bi,bj)*maskInC(i,j-1,bi,bj) #endif ENDDO ENDDO ENDIF ENDDO ENDDO IF ( calcMeanResid ) THEN _GLOBAL_SUM_RL( residU, myThid ) _GLOBAL_SUM_RL( residV, myThid ) C scale residuals by globalArea so that they do not get ridiculously large IF ( residU.GT.0. ) residU = SQRT(residU/globalArea) IF ( residV.GT.0. ) residV = SQRT(residV/globalArea) ENDIF C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| #endif /* SEAICE_ALLOW_DYNAMICS */ #endif /* SEAICE_CGRID */ RETURN END CBOP C !ROUTINE: SEAICE_LSR_CALC_COEFFS C !INTERFACE: SUBROUTINE SEAICE_LSR_CALC_COEFFS( I etaPlusZeta, zetaMinusEta, etaZloc, zetaZloc, dragSym, O AU, BU, CU, AV, BV, CV, uRt1, uRt2, vRt1, vRt2, I iMin, iMax, jMin, jMax, myTime, myIter, myThid ) C !DESCRIPTION: \bv C *==========================================================* C | S/R SEAICE_LSR_CALC_COEFFS C | o Calculate coefficient matrix for LSR solver C *==========================================================* C \ev C !USES: IMPLICIT NONE C === Global variables === #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" #include "SEAICE_SIZE.h" #include "SEAICE_PARAMS.h" #include "SEAICE.h" C !INPUT/OUTPUT PARAMETERS: C === Routine arguments === C myTime :: Simulation time C myIter :: Simulation timestep number C myThid :: my Thread Id. number _RL myTime INTEGER myIter INTEGER myThid INTEGER iMin, iMax, jMin, jMax C combinations of bulk and shear viscosity at C points C abbreviations _RL etaPlusZeta (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL zetaMinusEta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C shear and bulk viscosity at Z points _RL etaZloc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL zetaZloc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C symmetric ice-ocean stress coefficient _RL dragSym (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C coefficients of ice velocities in coefficient matrix C for both U and V-equation C diagonals of coefficient matrices _RL AU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL BU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL CU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL AV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL BV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL CV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C coefficients for lateral points, u(j+/-1) _RL uRt1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL uRt2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C coefficients for lateral points, v(i+/-1) _RL vRt1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL vRt2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) CEOP #ifdef SEAICE_CGRID #ifdef SEAICE_ALLOW_DYNAMICS C !LOCAL VARIABLES: C === Local variables === C i,j,bi,bj :: Loop counters INTEGER i, j, bi, bj _RL hFacM, hFacP C backward difference extrapolation factor (includes 1/deltaT) _RL bdfAlphaOverDt C strongly implicit coupling factor _RL strImpCplFac C fractional area at velocity points _RL areaW(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL areaS(1-OLx:sNx+OLx,1-OLy:sNy+OLy) C coefficients of ice velocities in coefficient matrix C for both U and V-equation C XX: double derivative in X C YY: double derivative in Y C XM: metric term with derivative in X C YM: metric term with derivative in Y _RL UXX (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL UYY (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL UXM (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL UYM (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL VXX (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL VYY (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL VXM (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL VYM (1-OLx:sNx+OLx,1-OLy:sNy+OLy) C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| C backward difference extrapolation factor bdfAlphaOverDt = 1. _d 0 IF ( SEAICEuseBDF2 ) THEN IF ( myIter.EQ.nIter0 .AND. SEAICEmomStartBDF.EQ.0 ) THEN bdfAlphaOverDt = 1. _d 0 ELSE bdfAlphaOverDt = 1.5 _d 0 ENDIF ENDIF C bdfAlphaOverDt = bdfAlphaOverDt / SEAICE_deltaTdyn C strImpCplFac = 0. _d 0 IF ( SEAICEuseStrImpCpl ) strImpCplFac = 1. _d 0 DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) IF ( SEAICEscaleSurfStress ) THEN DO J=jMin,jMax DO I=iMin,iMax areaW(I,J) = 0.5 _d 0*(AREA(I,J,bi,bj)+AREA(I-1,J,bi,bj)) areaS(I,J) = 0.5 _d 0*(AREA(I,J,bi,bj)+AREA(I,J-1,bi,bj)) ENDDO ENDDO ELSE DO J=jMin,jMax DO I=iMin,iMax areaW(I,J) = 1. _d 0 areaS(I,J) = 1. _d 0 ENDDO ENDDO ENDIF C C coefficients of uIce(I,J) and vIce(I,J) belonging to ... C DO J=jMin,jMax DO I=iMin-1,iMax C ... d/dx (eta+zeta) d/dx u UXX(I,J) = _dyF(I,J,bi,bj) * etaPlusZeta(I,J,bi,bj) & * _recip_dxF(I,J,bi,bj) C ... d/dx (zeta-eta) k1 u UXM(I,J) = _dyF(I,J,bi,bj) * zetaMinusEta(I,J,bi,bj) & * k1AtC(I,J,bi,bj) * 0.5 _d 0 ENDDO ENDDO DO J=jMin,jMax+1 DO I=iMin,iMax C ... d/dy eta d/dy u UYY(I,J) = _dxV(I,J,bi,bj) * & ( etaZloc(I,J,bi,bj) + strImpCplFac*zetaZloc(I,J,bi,bj) ) & * _recip_dyU(I,J,bi,bj) C ... d/dy eta k2 u UYM(I,J) = _dxV(I,J,bi,bj) * etaZloc(I,J,bi,bj) & * k2AtZ(I,J,bi,bj) * 0.5 _d 0 ENDDO ENDDO DO J=jMin,jMax DO I=iMin,iMax+1 C ... d/dx eta dv/dx VXX(I,J) = _dyU(I,J,bi,bj) * & ( etaZloc(I,J,bi,bj) + strImpCplFac*zetaZloc(I,J,bi,bj) ) & * _recip_dxV(I,J,bi,bj) C ... d/dx eta k1 v VXM(I,J) = _dyU(I,J,bi,bj) * etaZloc(I,J,bi,bj) & * k1AtZ(I,J,bi,bj) * 0.5 _d 0 ENDDO ENDDO DO J=jMin-1,jMax DO I=iMin,iMax C ... d/dy eta+zeta dv/dy VYY(I,J) = _dxF(I,J,bi,bj) * etaPlusZeta(I,J,bi,bj) & * _recip_dyF(I,J,bi,bj) C ... d/dy (zeta-eta) k2 v VYM(I,J) = _dxF(I,J,bi,bj) * zetaMinusEta(I,J,bi,bj) & * k2AtC(I,J,bi,bj) * 0.5 _d 0 ENDDO ENDDO C-- Coefficients for solving uIce : C assemble coefficient matrix, beware of sign convention: because this C is the left hand side we calculate -grad(sigma), but the coefficients C of U(I,J+/-1) are counted on the right hand side DO J=jMin,jMax DO I=iMin,iMax C coefficients for UICE(I-1,J) AU(I,J,bi,bj)= ( - UXX(I-1,J) + UXM(I-1,J) ) & * seaiceMaskU(I,J,bi,bj) C coefficients for UICE(I+1,J) CU(I,J,bi,bj)= ( - UXX(I ,J) - UXM(I ,J) ) & * seaiceMaskU(I,J,bi,bj) C coefficients for UICE(I,J) BU(I,J,bi,bj)=(ONE - seaiceMaskU(I,J,bi,bj)) + & ( UXX(I-1,J) + UXX(I,J) + UYY(I,J+1) + UYY(I,J) & + UXM(I-1,J) - UXM(I,J) + UYM(I,J+1) - UYM(I,J) & ) * seaiceMaskU(I,J,bi,bj) C coefficients of uIce(I,J-1) uRt1(I,J,bi,bj)= UYY(I,J ) + UYM(I,J ) C coefficients of uIce(I,J+1) uRt2(I,J,bi,bj)= UYY(I,J+1) - UYM(I,J+1) ENDDO ENDDO C apply boundary conditions according to slip factor C for no slip, set u on boundary to zero: u(j+/-1)=-u(j) C for the free slip case sigma_12 = 0 DO J=jMin,jMax DO I=iMin,iMax hFacM = seaiceMaskU(I,J-1,bi,bj) hFacP = seaiceMaskU(I,J+1,bi,bj) C copy contributions to coefficient of U(I,J) C beware of sign convection: uRt1/2 have the opposite sign convention C than BU, hence the minus sign BU(I,J,bi,bj)=BU(I,J,bi,bj) + seaiceMaskU(I,J,bi,bj) * & ( ( 1. _d 0 - hFacM ) * ( UYY(I ,J ) + UYM(I ,J ) ) & + ( 1. _d 0 - hFacP ) * ( UYY(I ,J+1) - UYM(I ,J+1) ) ) C reset coefficients of U(I,J-1) and U(I,J+1) uRt1(I,J,bi,bj) = uRt1(I,J,bi,bj) * hFacM uRt2(I,J,bi,bj) = uRt2(I,J,bi,bj) * hFacP ENDDO ENDDO C now we need to normalize everything by the grid cell area DO J=jMin,jMax DO I=iMin,iMax AU(I,J,bi,bj) = AU(I,J,bi,bj) * recip_rAw(I,J,bi,bj) CU(I,J,bi,bj) = CU(I,J,bi,bj) * recip_rAw(I,J,bi,bj) C here we need to add the contribution from the time derivative (in C bdfAlphaOverDt) and the symmetric drag term; must be done after C normalizing BU(I,J,bi,bj) = BU(I,J,bi,bj) * recip_rAw(I,J,bi,bj) & + seaiceMaskU(I,J,bi,bj) * & ( bdfAlphaOverDt*seaiceMassU(I,J,bi,bj) & + 0.5 _d 0 * ( dragSym(I, J,bi,bj) & + dragSym(I-1,J,bi,bj) )*areaW(I,J) & ) uRt1(I,J,bi,bj) = uRt1(I,J,bi,bj) * recip_rAw(I,J,bi,bj) uRt2(I,J,bi,bj) = uRt2(I,J,bi,bj) * recip_rAw(I,J,bi,bj) ENDDO ENDDO C-- Coefficients for solving uIce : C assemble coefficient matrix, beware of sign convention: because this C is the left hand side we calculate -grad(sigma), but the coefficients C of V(I+/-1,J) are counted on the right hand side DO J=jMin,jMax DO I=iMin,iMax C coefficients for VICE(I,J-1) AV(I,J,bi,bj)=( - VYY(I,J-1) + VYM(I,J-1) & ) * seaiceMaskV(I,J,bi,bj) C coefficients for VICE(I,J+1) CV(I,J,bi,bj)=( - VYY(I,J ) - VYM(I,J ) & ) * seaiceMaskV(I,J,bi,bj) C coefficients for VICE(I,J) BV(I,J,bi,bj)= (ONE - seaiceMaskV(I,J,bi,bj)) + & ( VXX(I,J) + VXX(I+1,J) + VYY(I,J) + VYY(I,J-1) & - VXM(I,J) + VXM(I+1,J) - VYM(I,J) + VYM(I,J-1) & ) * seaiceMaskV(I,J,bi,bj) C coefficients for V(I-1,J) vRt1(I,J,bi,bj) = VXX(I ,J) + VXM(I ,J) C coefficients for V(I+1,J) vRt2(I,J,bi,bj) = VXX(I+1,J) - VXM(I+1,J) ENDDO ENDDO C apply boundary conditions according to slip factor C for no slip, set u on boundary to zero: v(i+/-1)=-v(i) C for the free slip case sigma_12 = 0 DO J=jMin,jMax DO I=iMin,iMax hFacM = seaiceMaskV(i-1,j,bi,bj) hFacP = seaiceMaskV(i+1,j,bi,bj) C copy contributions to coefficient of V(I,J) C beware of sign convection: vRt1/2 have the opposite sign convention C than BV, hence the minus sign BV(I,J,bi,bj)=BV(I,J,bi,bj) + seaiceMaskV(I,J,bi,bj) * & ( ( 1. _d 0 - hFacM ) * ( VXX(I ,J) + VXM(I ,J) ) & + ( 1. _d 0 - hFacP ) * ( VXX(I+1,J) - VXM(I+1,J) ) ) C reset coefficients of V(I-1,J) and V(I+1,J) vRt1(I,J,bi,bj) = vRt1(I,J,bi,bj) * hFacM vRt2(I,J,bi,bj) = vRt2(I,J,bi,bj) * hFacP ENDDO ENDDO C now we need to normalize everything by the grid cell area DO J=jMin,jMax DO I=iMin,iMax AV(I,J,bi,bj) = AV(I,J,bi,bj) * recip_rAs(I,J,bi,bj) CV(I,J,bi,bj) = CV(I,J,bi,bj) * recip_rAs(I,J,bi,bj) C here we need add the contribution from the time derivative (in C bdfAlphaOverDt) and the symmetric drag term; must be done after C normalizing BV(I,J,bi,bj) = BV(I,J,bi,bj) * recip_rAs(I,J,bi,bj) & + seaiceMaskV(I,J,bi,bj) * & ( bdfAlphaOverDt*seaiceMassV(I,J,bi,bj) & + 0.5 _d 0 * ( dragSym(I,J, bi,bj) & + dragSym(I,J-1,bi,bj) )*areaS(I,J) & ) vRt1(I,J,bi,bj) = vRt1(I,J,bi,bj) * recip_rAs(I,J,bi,bj) vRt2(I,J,bi,bj) = vRt2(I,J,bi,bj) * recip_rAs(I,J,bi,bj) ENDDO ENDDO IF ( ( useCubedSphereExchange .AND. & (SEAICE_OLx.GT.0 .OR. SEAICE_OLy.GT.0) ) .OR. & SEAICEscaleSurfStress ) THEN C this is clearly a hack: make sure that diagonals BU/BV are non-zero. C In my experience we only need to catch the case of SEAICE_OLx/y=OLx/y-2, C but for safety lets call this routine whenever SEAICE_OLx/y > 0 C & (SEAICE_OLx.EQ.OLx-2 .OR. SEAICE_OLy.EQ.OLy-2) ) THEN C When scaling the surface ice-ocean stress by AREA, then there will C be many cases of zero diagonal entries. DO J=jMin,jMax DO I=iMin,iMax IF (BU(I,J,bi,bj).EQ.0 _d 0) BU(I,J,bi,bj) = 1. _d 0 IF (BV(I,J,bi,bj).EQ.0 _d 0) BV(I,J,bi,bj) = 1. _d 0 ENDDO ENDDO ENDIF C bi/bj-loops ENDDO ENDDO RETURN END C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP C !ROUTINE: SEAICE_LSR_RHSU C !INTERFACE: SUBROUTINE SEAICE_LSR_RHSU( I zetaMinusEta, etaPlusZeta, etaZloc, zetaZloc, pressLoc, I uIceC, vIceC, U rhsU, I iMin, iMax, jMin, jMax, bi, bj, myThid ) C !DESCRIPTION: \bv C *==========================================================* C | S/R SEAICE_LSR_RHSU C | o Set up stress contribution to rhs of u-momentum eq. C *==========================================================* C \ev C !USES: IMPLICIT NONE #include "SIZE.h" #include "EEPARAMS.h" #include "GRID.h" #include "SEAICE_SIZE.h" #include "SEAICE_PARAMS.h" #include "SEAICE.h" C !INPUT/OUTPUT PARAMETERS: C === Routine arguments === C myThid :: my Thread Id. number INTEGER myThid INTEGER iMin, iMax, jMin, jMax, bi, bj C _RL zetaMinusEta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL etaPlusZeta (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL etaZloc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL zetaZloc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL pressloc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL uIceC (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL vIceC (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL rhsU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C !LOCAL VARIABLES: C === Local variables === C i,j :: Loop counters INTEGER I,J INTEGER kSrf _RL hFacM _RL sig11(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL sig12(1-OLx:sNx+OLx,1-OLy:sNy+OLy) CEOP kSrf = 1 DO J=1-OLy,sNy+OLy DO I=1-OLx,sNx+OLx sig11(I,J) = 0. _d 0 sig12(I,J) = 0. _d 0 ENDDO ENDDO C contribution of sigma11 to rhs DO J=jMin,jMax DO I=iMin-1,iMax sig11(I,J) = zetaMinusEta(I,J,bi,bj) & * ( vIceC(I,J+1,bi,bj) - vIceC(I,J,bi,bj) ) & * _recip_dyF(I,J,bi,bj) & + etaPlusZeta(I,J,bi,bj) * k2AtC(I,J,bi,bj) & * 0.5 _d 0 * ( vIceC(I,J+1,bi,bj) + vIceC(I,J,bi,bj) ) & - 0.5 _d 0 * pressLoc(I,J,bi,bj) ENDDO ENDDO C contribution of sigma12 to rhs of u-equation DO J=jMin,jMax+1 DO I=iMin,iMax hFacM = seaiceMaskV(I,J,bi,bj) - seaiceMaskV(I-1,J,bi,bj) sig12(I,J) = etaZloc(I,J,bi,bj) * ( & ( vIceC(I,J,bi,bj) - vIceC(I-1,J,bi,bj) ) & * _recip_dxV(I,J,bi,bj) & - k1AtZ(I,J,bi,bj) & * 0.5 _d 0 * ( vIceC(I,J,bi,bj) + vIceC(I-1,J,bi,bj) ) & ) C free slip conditions (sig12=0) are taken care of by masking sig12 & *maskC(I ,J ,kSrf,bi,bj)*maskC(I-1,J ,kSrf,bi,bj) & *maskC(I ,J-1,kSrf,bi,bj)*maskC(I-1,J-1,kSrf,bi,bj) C no slip boundary conditions (v(i-1)=-v(i)) C v(i)+v(i-1) = 0 is also taken care of by masking sig12, so that we C only need to deal with v(i)-v(i-1) & + etaZloc(I,J,bi,bj) * _recip_dxV(I,J,bi,bj) & * ( vIceC(I,J,bi,bj) + vIceC(I-1,J,bi,bj) ) & * hFacM * 2. _d 0 ENDDO ENDDO IF ( SEAICEuseStrImpCpl ) THEN C strictly speaking, this is not a contribution to sig12, but for C convenience we add the explicit term -zetaZ*du/dy here; C we do not have to add any metric terms, because they cancel. DO J=jMin,jMax+1 DO I=iMin,iMax hFacM = seaiceMaskV(I,J,bi,bj) - seaiceMaskV(I-1,J,bi,bj) sig12(I,J) = sig12(I,J) - zetaZloc(I,J,bi,bj) * ( & ( uIceC(I,J,bi,bj) - uIceC(I,J-1,bi,bj) ) & * _recip_dyU(I,J,bi,bj) & ) C free slip conditions (sig12=0) are taken care of by masking sig12 & *maskC(I ,J ,kSrf,bi,bj)*maskC(I-1,J ,kSrf,bi,bj) & *maskC(I ,J-1,kSrf,bi,bj)*maskC(I-1,J-1,kSrf,bi,bj) C no slip boundary conditions (u(j-1)=-u(j)) C u(j)+u(j-1) = 0 is also taken care of by masking sig12, so that we C only need to deal with u(j)-u(j-1) & - zetaZloc(I,J,bi,bj) * _recip_dyU(I,J,bi,bj) & * ( uIceC(I,J,bi,bj) + uIceC(I,J-1,bi,bj) ) & * hFacM * 2. _d 0 ENDDO ENDDO ENDIF DO J=jMin,jMax DO I=iMin,iMax C contribution to the rhs part of grad(sigma)_x rhsU(I,J,bi,bj) = rhsU(I,J,bi,bj) & + recip_rAw(I,J,bi,bj) * seaiceMaskU(I,J,bi,bj) * & ( _dyF(I ,J ,bi,bj)*sig11(I ,J ) & - _dyF(I-1,J ,bi,bj)*sig11(I-1,J ) & + _dxV(I ,J+1,bi,bj)*sig12(I ,J+1) & - _dxV(I ,J ,bi,bj)*sig12(I ,J ) ) ENDDO ENDDO RETURN END C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP C !ROUTINE: SEAICE_LSR_RHSV C !INTERFACE: SUBROUTINE SEAICE_LSR_RHSV( I zetaMinusEta, etaPlusZeta, etaZloc, zetaZloc, pressLoc, I uIceC, vIceC, U rhsV, I iMin, iMax, jMin, jMax, bi, bj, myThid ) C !DESCRIPTION: \bv C *==========================================================* C | S/R SEAICE_LSR_RHSV C | o Set up stress contribution to rhs of v-momentum eq. C *==========================================================* C \ev C !USES: IMPLICIT NONE #include "SIZE.h" #include "EEPARAMS.h" #include "GRID.h" #include "SEAICE_SIZE.h" #include "SEAICE_PARAMS.h" #include "SEAICE.h" C !INPUT/OUTPUT PARAMETERS: C === Routine arguments === C myThid :: my Thread Id. number INTEGER myThid INTEGER iMin, iMax, jMin, jMax, bi, bj C _RL zetaMinusEta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL etaPlusZeta (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL etaZloc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL zetaZloc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL pressloc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL uIceC (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL vIceC (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL rhsV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C !LOCAL VARIABLES: C === Local variables === C i,j :: Loop counters INTEGER I,J INTEGER kSrf _RL hFacM _RL sig22(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL sig12(1-OLx:sNx+OLx,1-OLy:sNy+OLy) CEOP kSrf = 1 DO J = 1-Oly,sNy+Oly DO I = 1-Olx,sNx+Olx sig22(I,J) = 0. _d 0 sig12(I,J) = 0. _d 0 ENDDO ENDDO C contribution of sigma22 to rhs DO J=jMin-1,jMax DO I=iMin,iMax sig22(I,J) = zetaMinusEta(I,J,bi,bj) & * ( uIceC(I+1,J,bi,bj) - uIceC(I,J,bi,bj) ) & * _recip_dxF(I,J,bi,bj) & + etaPlusZeta(I,J,bi,bj) * k1AtC(I,J,bi,bj) & * 0.5 _d 0 * ( uIceC(I+1,J,bi,bj) + uIceC(I,J,bi,bj) ) & - 0.5 _d 0 * pressLoc(I,J,bi,bj) ENDDO ENDDO C contribution of sigma12 to rhs of v-equation DO J=jMin,jMax DO I=iMin,iMax+1 hFacM = seaiceMaskU(i,j,bi,bj) - seaiceMaskU(i,j-1,bi,bj) sig12(I,J) = etaZloc(I,J,bi,bj) * ( & ( uIceC(I,J,bi,bj) - uIceC(I,J-1,bi,bj) ) & * _recip_dyU(I,J,bi,bj) & - k2AtZ(I,J,bi,bj) & * 0.5 _d 0 * ( uIceC(I,J,bi,bj) + uIceC(I,J-1,bi,bj) ) & ) C free slip conditions (sig12=0) are taken care of by masking sig12, & *maskC(I ,J ,kSrf,bi,bj)*maskC(I-1,J ,kSrf,bi,bj) & *maskC(I ,J-1,kSrf,bi,bj)*maskC(I-1,J-1,kSrf,bi,bj) C no slip boundary conditions (u(j-1)=-u(j)) C u(j)+u(j-1) = 0 is also taken care of by masking sig12, so that we C only need to deal with u(j)-u(j-1) & + etaZloc(I,J,bi,bj) * _recip_dyU(I,J,bi,bj) & * ( uIceC(I,J,bi,bj) + uIceC(I,J-1,bi,bj) ) & * hFacM * 2. _d 0 ENDDO ENDDO IF ( SEAICEuseStrImpCpl ) THEN C strictly speaking, this is not a contribution to sig12, but for C convenience we add the explicit term -zetaZ*dv/dx here; C we do not have to add any metric terms, because they cancel. DO J=jMin,jMax DO I=iMin,iMax+1 hFacM = seaiceMaskU(i,j,bi,bj) - seaiceMaskU(i,j-1,bi,bj) sig12(I,J) = sig12(I,J) - zetaZloc(I,J,bi,bj) * ( & ( vIceC(I,J,bi,bj) - vIceC(I-1,J,bi,bj) ) & * _recip_dxV(I,J,bi,bj) & ) C free slip conditions (sig12=0) are taken care of by masking sig12 & *maskC(I ,J ,kSrf,bi,bj)*maskC(I-1,J ,kSrf,bi,bj) & *maskC(I ,J-1,kSrf,bi,bj)*maskC(I-1,J-1,kSrf,bi,bj) C no slip boundary conditions (v(i-1)=-v(i)) C v(i)+v(i-1) = 0 is also taken care of by masking sig12, so that we C only need to deal with v(i)-v(i-1) & - zetaZloc(I,J,bi,bj) * _recip_dxV(I,J,bi,bj) & * ( vIceC(I,J,bi,bj) + vIceC(I-1,J,bi,bj) ) & * hFacM * 2. _d 0 ENDDO ENDDO ENDIF DO J=jMin,jMax DO I=iMin,iMax C contribution to the rhs part of grad(sigma)_y rhsV(I,J,bi,bj) = rhsV(I,J,bi,bj) & + recip_rAs(I,J,bi,bj) * seaiceMaskV(I,J,bi,bj) * & ( _dyU(I+1,J ,bi,bj) * sig12(I+1,J ) & - _dyU(I ,J ,bi,bj) * sig12(I ,J ) & + _dxF(I ,J ,bi,bj) * sig22(I ,J ) & - _dxF(I ,J-1,bi,bj) * sig22(I ,J-1) ) ENDDO ENDDO RETURN END C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP C !ROUTINE: SEAICE_LSR_TRIDIAGU C !INTERFACE: SUBROUTINE SEAICE_LSR_TRIDIAGU( I AU, BU, CU, uRt1, uRt2, rhsU, uTmp, seaiceMaskU, WFAU, U uIce, I iMin, iMax, jMin, jMax, bi, bj, myTime, myIter, myThid ) C !DESCRIPTION: \bv C *==========================================================* C | S/R SEAICE_LSR_TRIDIAGU C | o Solve tridiagonal problem in uIce for LSR solver C *==========================================================* C \ev C !USES: IMPLICIT NONE #include "SIZE.h" C !INPUT/OUTPUT PARAMETERS: C === Routine arguments === C myTime :: Simulation time C myIter :: Simulation timestep number C myThid :: my Thread Id. number _RL myTime INTEGER myIter INTEGER myThid INTEGER iMin, iMax, jMin, jMax, bi, bj C diagonals of coefficient matrices _RL AU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL BU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL CU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C coefficients for lateral points, v(j+/-1) _RL uRt1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL uRt2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C right hand side _RL rhsU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C velocity of previous iteration step _RL uTmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C on entry: first guess and on exit: solution _RL uIce(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C mask _RL seaiceMaskU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C relaxation factor _RL WFAU CEOP C !LOCAL VARIABLES: C === Local variables === C i,j :: Loop counters INTEGER I,J,IM INTEGER jMinLoc, jStep #ifdef SEAICE_LSR_ZEBRA INTEGER K PARAMETER ( jStep = 2 ) #else PARAMETER ( jStep = 1 ) #endif /* SEAICE_LSR_ZEBRA */ _RL AA3, bet _RL URT(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL CUU(1-OLx:sNx+OLx,1-OLy:sNy+OLy) C Need to initialize local arrays DO J = 1-OLy,sNy+OLy DO I = 1-OLx,sNx+OLx URT(I,J) = 0. _d 0 CUU(I,J) = 0. _d 0 ENDDO ENDDO #ifdef SEAICE_LSR_ZEBRA DO K=0,1 jMinLoc = jMin+K #else jMinLoc = jMin #endif /* SEAICE_LSR_ZEBRA */ DO J=jMinLoc,jMax,jStep DO I=iMin,iMax AA3 = 0. _d 0 IF (I.EQ.iMin) AA3 = AA3 - AU(I,J,bi,bj)*uIce(I-1,J,bi,bj) IF (I.EQ.iMax) AA3 = AA3 - CU(I,J,bi,bj)*uIce(I+1,J,bi,bj) URT(I,J)=rhsU(I,J,bi,bj) & + AA3 & + uRt1(I,J,bi,bj)*uIce(I,J-1,bi,bj) & + uRt2(I,J,bi,bj)*uIce(I,J+1,bi,bj) URT(I,J)=URT(I,J)* seaiceMaskU(I,J,bi,bj) ENDDO C beginning of tridiagnoal solver DO I=iMin,iMax CUU(I,J)=CU(I,J,bi,bj) ENDDO CUU(iMin,J)=CUU(iMin,J)/BU(iMin,J,bi,bj) URT(iMin,J)=URT(iMin,J)/BU(iMin,J,bi,bj) #ifdef SEAICE_VECTORIZE_LSR ENDDO C start a new loop with reversed order to support automatic vectorization CADJ loop = sequential DO I=iMin+1,iMax IM=I-1 DO J=jMinLoc,jMax,jStep #else /* do not SEAICE_VECTORIZE_LSR */ CADJ loop = sequential DO I=iMin+1,iMax IM=I-1 #endif /* SEAICE_VECTORIZE_LSR */ bet = BU(I,J,bi,bj)-AU(I,J,bi,bj)*CUU(IM,J) CUU(I,J) = CUU(I,J)/bet URT(I,J) = (URT(I,J)-AU(I,J,bi,bj)*URT(IM,J))/bet ENDDO #ifdef SEAICE_VECTORIZE_LSR ENDDO CADJ loop = sequential DO I=iMin,iMax-1 IM=sNx-I DO J=jMinLoc,jMax,jStep #else CADJ loop = sequential DO I=iMin,iMax-1 IM=sNx-I #endif /* SEAICE_VECTORIZE_LSR */ URT(IM,J)=URT(IM,J)-CUU(IM,J)*URT(IM+1,J) ENDDO #ifdef SEAICE_VECTORIZE_LSR ENDDO #endif /* SEAICE_VECTORIZE_LSR */ C end of tridiagnoal solver, solution is URT C relaxation with WFAU #ifdef SEAICE_VECTORIZE_LSR C go back to original order DO J=jMinLoc,jMax,jStep #endif /* SEAICE_VECTORIZE_LSR */ DO I=iMin,iMax uIce(I,J,bi,bj)=uTmp(I,J,bi,bj) & +WFAU*(URT(I,J)-uTmp(I,J,bi,bj)) ENDDO ENDDO #ifdef SEAICE_LSR_ZEBRA ENDDO #endif /* SEAICE_LSR_ZEBRA */ RETURN END C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP C !ROUTINE: SEAICE_LSR_TRIDIAGV C !INTERFACE: SUBROUTINE SEAICE_LSR_TRIDIAGV( I AV, BV, CV, vRt1, vRt2, rhsV, vTmp, seaiceMaskV, WFAV, U vIce, I iMin, iMax, jMin, jMax, bi, bj, myTime, myIter, myThid ) C !DESCRIPTION: \bv C *==========================================================* C | S/R SEAICE_LSR_TRIDIAGV C | o Solve tridiagonal problem in vIce for LSR solver C *==========================================================* C \ev C !USES: IMPLICIT NONE #include "SIZE.h" C !INPUT/OUTPUT PARAMETERS: C === Routine arguments === C myTime :: Simulation time C myIter :: Simulation timestep number C myThid :: my Thread Id. number _RL myTime INTEGER myIter INTEGER myThid INTEGER iMin, iMax, jMin, jMax, bi, bj C diagonals of coefficient matrices _RL AV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL BV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL CV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C coefficients for lateral points, v(j+/-1) _RL vRt1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL vRt2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C right hand side _RL rhsV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C velocity of previous iteration step _RL vTmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C on entry: first guess and on exit: solution _RL vIce(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C mask _RL seaiceMaskV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C relaxation factor _RL WFAV CEOP C !LOCAL VARIABLES: C === Local variables === C i,j :: Loop counters INTEGER I,J,JM INTEGER iMinLoc, iStep #ifdef SEAICE_LSR_ZEBRA INTEGER K PARAMETER ( iStep = 2 ) #else PARAMETER ( iStep = 1 ) #endif /* SEAICE_LSR_ZEBRA */ _RL AA3, bet _RL VRT(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL CVV(1-OLx:sNx+OLx,1-OLy:sNy+OLy) C Need to initialize local arrays DO J = 1-OLy,sNy+OLy DO I = 1-OLx,sNx+OLx VRT(I,J) = 0. _d 0 CVV(I,J) = 0. _d 0 ENDDO ENDDO #ifdef SEAICE_LSR_ZEBRA DO K=0,1 iMinLoc = iMin+K #else iMinLoc = iMin #endif /* SEAICE_LSR_ZEBRA */ #ifdef SEAICE_VECTORIZE_LSR DO J=jMin,jMax DO I=iMinLoc,iMax,iStep #else DO I=iMinLoc,iMax,iStep DO J=jMin,jMax #endif /* SEAICE_VECTORIZE_LSR */ AA3 = 0. _d 0 IF (J.EQ.jMin) AA3 = AA3 - AV(I,J,bi,bj)*vIce(I,J-1,bi,bj) IF (J.EQ.jMax) AA3 = AA3 - CV(I,J,bi,bj)*vIce(I,J+1,bi,bj) VRT(I,J)=rhsV(I,J,bi,bj) & + AA3 & + vRt1(I,J,bi,bj)*vIce(I-1,J,bi,bj) & + vRt2(I,J,bi,bj)*vIce(I+1,J,bi,bj) VRT(I,J)=VRT(I,J)* seaiceMaskV(I,J,bi,bj) C beginning of tridiagnoal solver CVV(I,J)=CV(I,J,bi,bj) ENDDO #ifdef SEAICE_VECTORIZE_LSR ENDDO DO I=iMinLoc,iMax,iStep #endif /* SEAICE_VECTORIZE_LSR */ CVV(I,jMin)=CVV(I,jMin)/BV(I,jMin,bi,bj) VRT(I,jMin)=VRT(I,jMin)/BV(I,jMin,bi,bj) #ifdef SEAICE_VECTORIZE_LSR ENDDO C start a new loop with reversed order to support automatic vectorization CADJ loop = sequential DO J=jMin+1,jMax JM=J-1 DO I=iMinLoc,iMax,iStep #else /* do not SEAICE_VECTORIZE_LSR */ CADJ loop = sequential DO J=jMin+1,jMax JM=J-1 #endif /* SEAICE_VECTORIZE_LSR */ bet = BV(I,J,bi,bj)-AV(I,J,bi,bj)*CVV(I,JM) CVV(I,J) = CVV(I,J)/bet VRT(I,J) = (VRT(I,J)-AV(I,J,bi,bj)*VRT(I,JM))/bet ENDDO #ifdef SEAICE_VECTORIZE_LSR ENDDO C go back to original order CADJ loop = sequential DO J=jMin,jMax-1 JM=sNy-J DO I=iMinLoc,iMax,iStep #else CADJ loop = sequential DO J=jMin,jMax-1 JM=sNy-J #endif /* SEAICE_VECTORIZE_LSR */ VRT(I,JM)=VRT(I,JM)-CVV(I,JM)*VRT(I,JM+1) ENDDO #ifdef SEAICE_VECTORIZE_LSR ENDDO C end of tridiagnoal solver, solution is VRT C relaxation with WFAV DO J=jMin,jMax DO I=iMinLoc,iMax,iStep #else DO J=jMin,jMax #endif /* SEAICE_VECTORIZE_LSR */ vIce(I,J,bi,bj)=vTmp(I,J,bi,bj) & +WFAV*(VRT(I,J)-vTmp(I,J,bi,bj)) ENDDO ENDDO #ifdef SEAICE_LSR_ZEBRA ENDDO #endif /* SEAICE_LSR_ZEBRA */ #endif /* SEAICE_ALLOW_DYNAMICS */ #endif /* SEAICE_CGRID */ RETURN END