C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_dynsolver.F,v 1.65 2017/05/09 02:54:16 jmc Exp $ C $Name: $ #include "SEAICE_OPTIONS.h" #ifdef ALLOW_AUTODIFF # include "AUTODIFF_OPTIONS.h" #endif CBOP C !ROUTINE: SEAICE_DYNSOLVER C !INTERFACE: SUBROUTINE SEAICE_DYNSOLVER( myTime, myIter, myThid ) C *==========================================================* C | SUBROUTINE SEAICE_DYNSOLVER C | o Ice dynamics using LSR solver C | Zhang and Hibler, JGR, 102, 8691-8702, 1997 C | or EVP explicit solver by Hunke and Dukowicz, JPO 27, C | 1849-1867 (1997) C *==========================================================* C | written by Martin Losch, March 2006 C *==========================================================* C \ev C !USES: IMPLICIT NONE C === Global variables === #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" #include "SURFACE.h" #include "DYNVARS.h" #include "FFIELDS.h" #include "SEAICE_SIZE.h" #include "SEAICE_PARAMS.h" #include "SEAICE.h" #ifdef ALLOW_AUTODIFF_TAMC # 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 C !FUNCTIONS: LOGICAL DIFFERENT_MULTIPLE EXTERNAL C !LOCAL VARIABLES: C === Local variables === C i,j,bi,bj :: Loop counters INTEGER i, j, bi, bj _RL phiSurf(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL mask_uice, mask_vice #ifdef ALLOW_AUTODIFF_TAMC _RL PSTAR #endif /* ALLOW_AUTODIFF_TAMC */ # ifdef ALLOW_AUTODIFF_TAMC C Following re-initialisation breaks some "artificial" AD dependencies C incured by IF (DIFFERENT_MULTIPLE ... statement PSTAR = SEAICE_strength DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1-OLy+1,sNy+OLy DO i=1-OLx+1,sNx+OLx PRESS0(i,j,bi,bj) = PSTAR*HEFF(i,j,bi,bj) & *EXP(-SEAICE_cStar*(ONE-AREA(i,j,bi,bj))) ZMAX(I,J,bi,bj) = SEAICE_zetaMaxFac*PRESS0(I,J,bi,bj) ZMIN(i,j,bi,bj) = SEAICE_zetaMin PRESS0(i,j,bi,bj) = PRESS0(i,j,bi,bj)*HEFFM(i,j,bi,bj) TAUX(i,j,bi,bj) = 0. _d 0 TAUY(i,j,bi,bj) = 0. _d 0 #ifdef SEAICE_ALLOW_FREEDRIFT uice_fd(i,j,bi,bj)= 0. _d 0 vice_fd(i,j,bi,bj)= 0. _d 0 #endif ENDDO ENDDO ENDDO ENDDO C CADJ STORE uice = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE vice = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE uicenm1 = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE vicenm1 = comlev1, key=ikey_dynamics, kind=isbyte # endif /* ALLOW_AUTODIFF_TAMC */ IF ( & DIFFERENT_MULTIPLE(SEAICE_deltaTdyn,myTime,SEAICE_deltaTtherm) & ) THEN # ifdef ALLOW_AUTODIFF_TAMC # ifdef SEAICE_ALLOW_EVP CADJ STORE press0 = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE zmax = comlev1, key=ikey_dynamics, kind=isbyte # endif # endif /* ALLOW_AUTODIFF_TAMC */ C-- NOW SET UP MASS PER UNIT AREA AND CORIOLIS TERM DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1-OLy+1,sNy+OLy DO i=1-OLx+1,sNx+OLx seaiceMassC(I,J,bi,bj)=SEAICE_rhoIce*HEFF(i,j,bi,bj) seaiceMassU(I,J,bi,bj)=SEAICE_rhoIce*HALF*( & HEFF(i,j,bi,bj) + HEFF(i-1,j ,bi,bj) ) seaiceMassV(I,J,bi,bj)=SEAICE_rhoIce*HALF*( & HEFF(i,j,bi,bj) + HEFF(i ,j-1,bi,bj) ) ENDDO ENDDO IF ( SEAICEaddSnowMass ) THEN DO j=1-OLy+1,sNy+OLy DO i=1-OLx+1,sNx+OLx seaiceMassC(I,J,bi,bj)=seaiceMassC(I,J,bi,bj) & + SEAICE_rhoSnow*HSNOW(i,j,bi,bj) seaiceMassU(I,J,bi,bj)=seaiceMassU(I,J,bi,bj) & + SEAICE_rhoSnow*HALF*( & HSNOW(i,j,bi,bj) + HSNOW(i-1,j ,bi,bj) ) seaiceMassV(I,J,bi,bj)=seaiceMassV(I,J,bi,bj) & + SEAICE_rhoSnow*HALF*( & HSNOW(i,j,bi,bj) + HSNOW(i ,j-1,bi,bj) ) ENDDO ENDDO ENDIF ENDDO ENDDO #ifndef ALLOW_AUTODIFF_TAMC IF ( SEAICE_maskRHS ) THEN C dynamic masking of areas with no ice, not recommended C and only kept for testing purposes DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1-OLy+1,sNy+OLy DO i=1-OLx+1,sNx+OLx seaiceMaskU(I,J,bi,bj)=AREA(i,j,bi,bj)+AREA(I-1,J,bi,bj) mask_uice=HEFFM(i,j,bi,bj)+HEFFM(i-1,j ,bi,bj) IF ( (seaiceMaskU(I,J,bi,bj) .GT. 0. _d 0) .AND. & (mask_uice .GT. 1.5 _d 0) ) THEN seaiceMaskU(I,J,bi,bj) = 1. _d 0 ELSE seaiceMaskU(I,J,bi,bj) = 0. _d 0 ENDIF seaiceMaskV(I,J,bi,bj)=AREA(i,j,bi,bj)+AREA(I,J-1,bi,bj) mask_vice=HEFFM(i,j,bi,bj)+HEFFM(i ,j-1,bi,bj) IF ( (seaiceMaskV(I,J,bi,bj) .GT. 0. _d 0) .AND. & (mask_vice .GT. 1.5 _d 0) ) THEN seaiceMaskV(I,J,bi,bj) = 1. _d 0 ELSE seaiceMaskV(I,J,bi,bj) = 0. _d 0 ENDIF ENDDO ENDDO ENDDO ENDDO CALL EXCH_UV_XY_RL( seaiceMaskU, seaiceMaskV, .FALSE., myThid ) ENDIF #endif /* ndef ALLOW_AUTODIFF_TAMC */ C-- NOW SET UP FORCING FIELDS C initialise fields DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx TAUX (I,J,bi,bj)= 0. _d 0 TAUY (I,J,bi,bj)= 0. _d 0 #ifdef ALLOW_AUTODIFF_TAMC # ifdef SEAICE_ALLOW_EVP stressDivergenceX(I,J,bi,bj) = 0. _d 0 stressDivergenceY(I,J,bi,bj) = 0. _d 0 # endif #endif ENDDO ENDDO ENDDO ENDDO # ifdef ALLOW_AUTODIFF_TAMC CADJ STORE uice = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE vice = comlev1, key=ikey_dynamics, kind=isbyte # endif /* ALLOW_AUTODIFF_TAMC */ C-- interface of dynamics with atmopheric forcing fields (wind/stress) CALL SEAICE_GET_DYNFORCING ( I uIce, vIce, AREA, O TAUX, TAUY, I myTime, myIter, myThid ) DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) C-- Compute surface pressure at z==0: C- use actual sea surface height for tilt computations DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx phiSurf(i,j) = Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj) ENDDO ENDDO #ifdef ATMOSPHERIC_LOADING C- add atmospheric loading and Sea-Ice loading IF ( useRealFreshWaterFlux ) THEN DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx phiSurf(i,j) = phiSurf(i,j) & + ( pload(i,j,bi,bj) & +sIceLoad(i,j,bi,bj)*gravity & )*recip_rhoConst ENDDO ENDDO ELSE DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx phiSurf(i,j) = phiSurf(i,j) & + pload(i,j,bi,bj)*recip_rhoConst ENDDO ENDDO ENDIF #endif /* ATMOSPHERIC_LOADING */ C-- basic forcing by wind stress IF ( SEAICEscaleSurfStress ) THEN DO j=1-OLy+1,sNy+OLy DO i=1-OLx+1,sNx+OLx FORCEX0(I,J,bi,bj)=TAUX(I,J,bi,bj) & * 0.5 _d 0*(AREA(I,J,bi,bj)+AREA(I-1,J,bi,bj)) FORCEY0(I,J,bi,bj)=TAUY(I,J,bi,bj) & * 0.5 _d 0*(AREA(I,J,bi,bj)+AREA(I,J-1,bi,bj)) ENDDO ENDDO ELSE DO j=1-OLy+1,sNy+OLy DO i=1-OLx+1,sNx+OLx FORCEX0(I,J,bi,bj)=TAUX(I,J,bi,bj) FORCEY0(I,J,bi,bj)=TAUY(I,J,bi,bj) ENDDO ENDDO ENDIF IF ( SEAICEuseTILT ) then DO j=1-OLy+1,sNy+OLy DO i=1-OLx+1,sNx+OLx C-- now add in tilt FORCEX0(I,J,bi,bj)=FORCEX0(I,J,bi,bj) & -seaiceMassU(I,J,bi,bj)*_recip_dxC(I,J,bi,bj) & *( phiSurf(i,j)-phiSurf(i-1,j) ) FORCEY0(I,J,bi,bj)=FORCEY0(I,J,bi,bj) & -seaiceMassV(I,J,bi,bj)* _recip_dyC(I,J,bi,bj) & *( phiSurf(i,j)-phiSurf(i,j-1) ) ENDDO ENDDO ENDIF CALL SEAICE_CALC_ICE_STRENGTH( bi, bj, myTime, myIter, myThid ) ENDDO ENDDO #ifdef SEAICE_ALLOW_DYNAMICS IF ( SEAICEuseDYNAMICS ) THEN #ifdef SEAICE_ALLOW_FREEDRIFT IF ( SEAICEuseFREEDRIFT .OR. SEAICEuseEVP & .OR. LSR_mixIniGuess.EQ.0 ) THEN CALL SEAICE_FREEDRIFT( myTime, myIter, myThid ) ENDIF IF ( SEAICEuseFREEDRIFT ) THEN 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_fd(i,j,bi,bj) vIce(i,j,bi,bj) = vIce_fd(i,j,bi,bj) stressDivergenceX(i,j,bi,bj) = 0. _d 0 stressDivergenceY(i,j,bi,bj) = 0. _d 0 ENDDO ENDDO ENDDO ENDDO ENDIF #endif /* SEAICE_ALLOW_FREEDRIFT */ #ifdef ALLOW_OBCS IF ( useOBCS ) THEN CALL OBCS_APPLY_UVICE( uIce, vIce, myThid ) ENDIF #endif /* ALLOW_OBCS */ #ifdef SEAICE_ALLOW_EVP # ifdef ALLOW_AUTODIFF_TAMC CADJ STORE uice = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE vice = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE uicenm1 = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE vicenm1 = comlev1, key=ikey_dynamics, kind=isbyte # endif /* ALLOW_AUTODIFF_TAMC */ IF ( SEAICEuseEVP ) THEN C Elastic-Viscous-Plastic solver, following Hunke (2001) CALL SEAICE_EVP( myTime, myIter, myThid ) ENDIF #endif /* SEAICE_ALLOW_EVP */ IF ( SEAICEuseLSR ) THEN C Picard solver with LSR scheme (Zhang-J/Hibler 1997), ported to a C-grid CALL SEAICE_LSR( myTime, myIter, myThid ) ENDIF #ifdef SEAICE_ALLOW_KRYLOV # ifdef ALLOW_AUTODIFF_TAMC STOP 'Adjoint does not work with Picard-Krylov solver.' else IF ( SEAICEuseKrylov ) THEN C Picard solver with Matrix-free Krylov solver (Lemieux et al. 2008) CALL SEAICE_KRYLOV( myTime, myIter, myThid ) ENDIF # endif /* ALLOW_AUTODIFF_TAMC */ #endif /* SEAICE_ALLOW_KRYLOV */ #ifdef SEAICE_ALLOW_JFNK # ifdef ALLOW_AUTODIFF_TAMC STOP 'Adjoint does not work with JFNK solver.' else IF ( SEAICEuseJFNK ) THEN C Jacobian-free Newton Krylov solver (Lemieux et al. 2010, 2012) CALL SEAICE_JFNK( myTime, myIter, myThid ) ENDIF # endif /* ALLOW_AUTODIFF_TAMC */ #endif /* SEAICE_ALLOW_JFNK */ C End of IF (SEAICEuseDYNAMICS ... ENDIF #endif /* SEAICE_ALLOW_DYNAMICS */ C End of IF (DIFFERENT_MULTIPLE ... ENDIF #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE uice = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE vice = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE stressDivergenceX = comlev1, CADJ & key=ikey_dynamics, kind=isbyte CADJ STORE stressDivergenceY = comlev1, CADJ & key=ikey_dynamics, kind=isbyte CADJ STORE DWATN = comlev1, key=ikey_dynamics, kind=isbyte #ifdef SEAICE_ALLOW_BOTTOMDRAG CADJ STORE CbotC = comlev1, key=ikey_dynamics, kind=isbyte #endif /* SEAICE_ALLOW_BOTTOMDRAG */ #endif /* ALLOW_AUTODIFF_TAMC */ C Calculate ocean surface stress CALL SEAICE_OCEAN_STRESS ( myTime, myIter, myThid ) #ifdef SEAICE_ALLOW_DYNAMICS #ifdef SEAICE_ALLOW_CLIPVELS IF ( SEAICEuseDYNAMICS .AND. SEAICE_clipVelocities) THEN #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE uice = comlev1, key=ikey_dynamics, kind=isbyte CADJ STORE vice = comlev1, key=ikey_dynamics, kind=isbyte #endif /* ALLOW_AUTODIFF_TAMC */ c Put a cap on ice velocity c limit velocity to 0.40 m s-1 to avoid potential CFL violations c in open water areas (drift of zero thickness ice) 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)= & MAX(MIN(uIce(i,j,bi,bj),0.40 _d +00),-0.40 _d +00) vIce(i,j,bi,bj)= & MAX(MIN(vIce(i,j,bi,bj),0.40 _d +00),-0.40 _d +00) ENDDO ENDDO ENDDO ENDDO ENDIF #endif /* SEAICE_ALLOW_CLIPVELS */ #endif /* SEAICE_ALLOW_DYNAMICS */ #endif /* SEAICE_CGRID */ RETURN END