pkg/seaice/dynsolver.F

809c36b928 Patr*0001 #include "SEAICE_OPTIONS.h"
8377b8ee87 Mart*0002 #ifdef ALLOW_EXF
                0003 # include "EXF_OPTIONS.h"
                0004 #endif
772b2ed80e Gael*0005 #ifdef ALLOW_AUTODIFF
                0006 # include "AUTODIFF_OPTIONS.h"
                0007 #endif
09510da3bb Dimi*0008
809c36b928 Patr*0009 CStartOfInterface
5d6433c292 Jean*0010       SUBROUTINE DYNSOLVER( myTime, myIter, myThid )
                0011 C     *==========================================================*
                0012 C     | SUBROUTINE DYNSOLVER                                     |
45315406aa Mart*0013 C     | o B-grid version of ice dynamics using LSR solver        |
                0014 C     |   Zhang and Hibler, JGR, 102, 8691-8702, 1997            |
5d6433c292 Jean*0015 C     *==========================================================*
                0016 C     *==========================================================*
809c36b928 Patr*0017       IMPLICIT NONE
09510da3bb Dimi*0018
809c36b928 Patr*0019 C     === Global variables ===
                0020 #include "SIZE.h"
                0021 #include "EEPARAMS.h"
                0022 #include "PARAMS.h"
4366d31d92 Mart*0023 #include "GRID.h"
55f4fb1c97 Mart*0024 #include "DYNVARS.h"
809c36b928 Patr*0025 #include "FFIELDS.h"
03c669d1ab Jean*0026 #include "SEAICE_SIZE.h"
809c36b928 Patr*0027 #include "SEAICE_PARAMS.h"
3f0f10fc37 Mart*0028 #include "SEAICE_GRID.h"
03c669d1ab Jean*0029 #include "SEAICE.h"
ae1fb66b64 Dimi*0030 #ifdef ALLOW_EXF
                0031 # include "EXF_FIELDS.h"
                0032 #endif
baa476eeba Dimi*0033 #ifdef ALLOW_AUTODIFF_TAMC
                0034 # include "tamc.h"
6060ec2938 Dimi*0035 #endif
baa476eeba Dimi*0036
809c36b928 Patr*0037 C     === Routine arguments ===
                0038 C     myTime - Simulation time
                0039 C     myIter - Simulation timestep number
                0040 C     myThid - Thread no. that called this routine.
                0041       _RL     myTime
                0042       INTEGER myIter
                0043       INTEGER myThid
                0044 CEndOfInterface
09510da3bb Dimi*0045
45315406aa Mart*0046 #ifdef SEAICE_BGRID_DYNAMICS
f3ce416a61 Jean*0047
                0048 #ifdef EXPLICIT_SSH_SLOPE
                0049 #include "SURFACE.h"
642c38b9c7 Jean*0050       _RL phiSurf(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
f3ce416a61 Jean*0051 #endif
                0052
809c36b928 Patr*0053 C     === Local variables ===
cee16b76ae Dimi*0054 C     i,j,bi,bj - Loop counters
809c36b928 Patr*0055
642c38b9c7 Jean*0056       INTEGER i, j, bi, bj
35fda33b05 Jean*0057       _RL RHOICE, RHOAIR
df93b38141 Mart*0058       _RL COSWIN
                0059       _RS SINWIN
642c38b9c7 Jean*0060       _RL ECCEN, ECM2, PSTAR, AAA
                0061       _RL U1, V1
809c36b928 Patr*0062
                0063       _RL COR_ICE    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,  nSx,nSy)
                0064
cee16b76ae Dimi*0065 C--   FIRST SET UP BASIC CONSTANTS
53092bcb42 Mart*0066       RHOICE=SEAICE_rhoIce
                0067       RHOAIR=SEAICE_rhoAir
                0068       ECCEN=SEAICE_eccen
cee16b76ae Dimi*0069       ECM2=ONE/(ECCEN**2)
                0070       PSTAR=SEAICE_strength
                0071
16f85413ea Mart*0072 C--   introduce turning angle (default is zero)
55f4fb1c97 Mart*0073       SINWIN=SIN(SEAICE_airTurnAngle*deg2rad)
                0074       COSWIN=COS(SEAICE_airTurnAngle*deg2rad)
                0075
                0076 C--   Compute proxy for geostrophic velocity,
                0077       DO bj=myByLo(myThid),myByHi(myThid)
                0078        DO bi=myBxLo(myThid),myBxHi(myThid)
                0079         DO j=0,sNy+1
                0080          DO i=0,sNx+1
8377b8ee87 Mart*0081           GWATX(i,j,bi,bj)=HALF*(uVel(i,j,KGEO(i,j,bi,bj),bi,bj)
                0082      &                         +uVel(i,j-1,KGEO(i,j,bi,bj),bi,bj))
                0083           GWATY(i,j,bi,bj)=HALF*(vVel(i,j,KGEO(i,j,bi,bj),bi,bj)
                0084      &                         +vVel(i-1,j,KGEO(i,j,bi,bj),bi,bj))
55f4fb1c97 Mart*0085 #ifdef SEAICE_DEBUG
                0086 c          write(*,'(2i4,2i2,f7.1,7f12.3)')
                0087 c     &     ,i,j,bi,bj,UVM(I,J,bi,bj)
                0088 c     &     ,GWATX(I,J,bi,bj),GWATY(I,J,bi,bj)
                0089 c     &     ,uVel(i+1,j,3,bi,bj),uVel(i+1,j+1,3,bi,bj)
                0090 c     &     ,vVel(i,j+1,3,bi,bj),vVel(i+1,j+1,3,bi,bj)
                0091 #endif
                0092          ENDDO
                0093         ENDDO
                0094        ENDDO
                0095       ENDDO
809c36b928 Patr*0096
cee16b76ae Dimi*0097 C--   NOW SET UP MASS PER UNIT AREA AND CORIOLIS TERM
                0098       DO bj=myByLo(myThid),myByHi(myThid)
                0099        DO bi=myBxLo(myThid),myBxHi(myThid)
0b5d39a040 Jean*0100         DO j=1-OLy,sNy+OLy
                0101          DO i=1-OLx,sNx+OLx
                0102           COR_ICE(i,j,bi,bj) = 0.
                0103          ENDDO
                0104         ENDDO
cee16b76ae Dimi*0105         DO j=1,sNy
                0106          DO i=1,sNx
8377b8ee87 Mart*0107           AMASS(i,j,bi,bj)=RHOICE*QUART*(
772590b63c Mart*0108      &          HEFF(i,j  ,bi,bj) + HEFF(i-1,j  ,bi,bj)
                0109      &         +HEFF(i,j-1,bi,bj) + HEFF(i-1,j-1,bi,bj) )
8377b8ee87 Mart*0110           COR_ICE(i,j,bi,bj)=AMASS(i,j,bi,bj) * _fCoriG(i,j,bi,bj)
cee16b76ae Dimi*0111          ENDDO
                0112         ENDDO
                0113        ENDDO
                0114       ENDDO
809c36b928 Patr*0115
cee16b76ae Dimi*0116 C--   NOW SET UP FORCING FIELDS
                0117
bda7c25f4a Dimi*0118 C--   Wind stress is computed on South-West B-grid U/V
                0119 C     locations from wind on tracer locations
                0120       DO bj=myByLo(myThid),myByHi(myThid)
                0121        DO bi=myBxLo(myThid),myBxHi(myThid)
                0122         DO j=1,sNy
                0123          DO i=1,sNx
8377b8ee87 Mart*0124           U1=QUART*(UWIND(i-1,j-1,bi,bj)+UWIND(i-1,j,bi,bj)
                0125      &             +UWIND(i  ,j-1,bi,bj)+UWIND(i  ,j,bi,bj))
                0126           V1=QUART*(VWIND(i-1,j-1,bi,bj)+VWIND(i-1,j,bi,bj)
                0127      &             +VWIND(i  ,j-1,bi,bj)+VWIND(i  ,j,bi,bj))
bda7c25f4a Dimi*0128           AAA=U1**2+V1**2
                0129           IF ( AAA .LE. SEAICE_EPS_SQ ) THEN
                0130              AAA=SEAICE_EPS
                0131           ELSE
                0132              AAA=SQRT(AAA)
                0133           ENDIF
09510da3bb Dimi*0134 C first ocean surface stress
8377b8ee87 Mart*0135           DAIRN(i,j,bi,bj)=RHOAIR*OCEAN_drag
bda7c25f4a Dimi*0136      &         *(2.70 _d 0+0.142 _d 0*AAA+0.0764 _d 0*AAA*AAA)
8377b8ee87 Mart*0137           WINDX(i,j,bi,bj)=DAIRN(i,j,bi,bj)*
                0138      &         (COSWIN*U1-SIGN(SINWIN, _fCori(i,j,bi,bj))*V1)
                0139           WINDY(i,j,bi,bj)=DAIRN(i,j,bi,bj)*
                0140      &         (SIGN(SINWIN, _fCori(i,j,bi,bj))*U1+COSWIN*V1)
09510da3bb Dimi*0141
                0142 C now ice surface stress
8377b8ee87 Mart*0143           IF ( YC(i,j,bi,bj) .LT. ZERO ) THEN
                0144            DAIRN(i,j,bi,bj) =
                0145      &          RHOAIR*(SEAICE_drag_south*AAA*AREA(i,j,bi,bj)
f834b21bef Dimi*0146      &          +OCEAN_drag*(2.70 _d 0+0.142 _d 0*AAA
8377b8ee87 Mart*0147      &          +0.0764 _d 0*AAA*AAA)*(ONE-AREA(i,j,bi,bj)))
f834b21bef Dimi*0148           ELSE
8377b8ee87 Mart*0149            DAIRN(i,j,bi,bj) =
                0150      &          RHOAIR*(SEAICE_drag*AAA*AREA(i,j,bi,bj)
f834b21bef Dimi*0151      &          +OCEAN_drag*(2.70 _d 0+0.142 _d 0*AAA
8377b8ee87 Mart*0152      &          +0.0764 _d 0*AAA*AAA)*(ONE-AREA(i,j,bi,bj)))
f834b21bef Dimi*0153           ENDIF
8377b8ee87 Mart*0154           FORCEX(i,j,bi,bj)=DAIRN(i,j,bi,bj)*
                0155      &         (COSWIN*U1-SIGN(SINWIN, _fCori(i,j,bi,bj))*V1)
                0156           FORCEY(i,j,bi,bj)=DAIRN(i,j,bi,bj)*
                0157      &         (SIGN(SINWIN, _fCori(i,j,bi,bj))*U1+COSWIN*V1)
cee16b76ae Dimi*0158          ENDDO
                0159         ENDDO
                0160        ENDDO
bda7c25f4a Dimi*0161       ENDDO
809c36b928 Patr*0162
                0163       DO bj=myByLo(myThid),myByHi(myThid)
                0164        DO bi=myBxLo(myThid),myBxHi(myThid)
f3ce416a61 Jean*0165 #ifdef EXPLICIT_SSH_SLOPE
                0166 C--   Compute surface pressure at z==0:
                0167 C-    use actual sea surface height for tilt computations
642c38b9c7 Jean*0168         DO j=1-OLy,sNy+OLy
                0169           DO i=1-OLx,sNx+OLx
f3ce416a61 Jean*0170             phiSurf(i,j) = Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
                0171           ENDDO
                0172         ENDDO
                0173 #ifdef ATMOSPHERIC_LOADING
                0174 C-    add atmospheric loading and Sea-Ice loading
                0175         IF ( useRealFreshWaterFlux ) THEN
642c38b9c7 Jean*0176           DO j=1-OLy,sNy+OLy
                0177            DO i=1-OLx,sNx+OLx
f3ce416a61 Jean*0178             phiSurf(i,j) = phiSurf(i,j)
                0179      &                   + ( pload(i,j,bi,bj)
0320e25227 Mart*0180      &                      +sIceLoad(i,j,bi,bj)*gravity*sIceLoadFac
f3ce416a61 Jean*0181      &                     )*recip_rhoConst
                0182            ENDDO
                0183           ENDDO
                0184         ELSE
642c38b9c7 Jean*0185           DO j=1-OLy,sNy+OLy
                0186            DO i=1-OLx,sNx+OLx
f3ce416a61 Jean*0187             phiSurf(i,j) = phiSurf(i,j)
                0188      &                   + pload(i,j,bi,bj)*recip_rhoConst
                0189            ENDDO
                0190           ENDDO
                0191         ENDIF
                0192 #endif /* ATMOSPHERIC_LOADING */
642c38b9c7 Jean*0193         DO j=1-OLy+1,sNy+OLy
                0194          DO i=1-OLx+1,sNx+OLx
f3ce416a61 Jean*0195 C--   NOW ADD IN TILT
8377b8ee87 Mart*0196           FORCEX(i,j,bi,bj)=FORCEX(i,j,bi,bj)
                0197      &      -AMASS(i,j,bi,bj)
ec0d7df165 Mart*0198      &         *( (phiSurf(i, j )-phiSurf(i-1, j ))*SIMaskU(i, j ,bi,bj)
                0199      &           +(phiSurf(i,j-1)-phiSurf(i-1,j-1))*SIMaskV(i,j-1,bi,bj)
8377b8ee87 Mart*0200      &          )*HALF*_recip_dxV(i,j,bi,bj)
                0201           FORCEY(i,j,bi,bj)=FORCEY(i,j,bi,bj)
                0202      &      -AMASS(i,j,bi,bj)
ec0d7df165 Mart*0203      &         *( (phiSurf( i ,j)-phiSurf( i ,j-1))*SIMaskV( i ,j,bi,bj)
                0204      &           +(phiSurf(i-1,j)-phiSurf(i-1,j-1))*SIMaskV(i-1,j,bi,bj)
8377b8ee87 Mart*0205      &          )*HALF*_recip_dyU(i,j,bi,bj)
f3ce416a61 Jean*0206 C NOW KEEP FORCEX0
8377b8ee87 Mart*0207           FORCEX0(i,j,bi,bj)=FORCEX(i,j,bi,bj)
                0208           FORCEY0(i,j,bi,bj)=FORCEY(i,j,bi,bj)
f3ce416a61 Jean*0209          ENDDO
                0210         ENDDO
                0211 #endif /* EXPLICIT_SSH_SLOPE */
642c38b9c7 Jean*0212         DO j=1-OLy,sNy+OLy
                0213          DO i=1-OLx,sNx+OLx
f3ce416a61 Jean*0214 #ifndef EXPLICIT_SSH_SLOPE
cee16b76ae Dimi*0215 C--   NOW ADD IN TILT
8377b8ee87 Mart*0216           FORCEX(i,j,bi,bj)=FORCEX(i,j,bi,bj)
                0217      &         -COR_ICE(i,j,bi,bj)*GWATY(i,j,bi,bj)
                0218           FORCEY(i,j,bi,bj)=FORCEY(i,j,bi,bj)
                0219      &         +COR_ICE(i,j,bi,bj)*GWATX(i,j,bi,bj)
09510da3bb Dimi*0220 C NOW KEEP FORCEX0
8377b8ee87 Mart*0221           FORCEX0(i,j,bi,bj)=FORCEX(i,j,bi,bj)
                0222           FORCEY0(i,j,bi,bj)=FORCEY(i,j,bi,bj)
f3ce416a61 Jean*0223 #endif /* EXPLICIT_SSH_SLOPE */
cee16b76ae Dimi*0224 C--   NOW SET UP ICE PRESSURE AND VISCOSITIES
8377b8ee87 Mart*0225           PRESS0(i,j,bi,bj)=PSTAR*HEFF(i,j,bi,bj)
ba6cfc5714 Mart*0226      &         *EXP(-SEAICE_cStar*(ONE-AREA(i,j,bi,bj)))
8e32c48b8f Mart*0227 CML          SEAICE_zMax(I,J,bi,bj)=(5.0 _d +12/(2.0 _d +04))*PRESS0(I,J,bi,bj)
                0228           SEAICE_zMax(i,j,bi,bj)=SEAICE_zetaMaxFac*PRESS0(i,j,bi,bj)
                0229 CML          SEAICE_zMin(I,J,bi,bj)=4.0 _d +08
                0230           SEAICE_zMin(i,j,bi,bj)=SEAICE_zetaMin
8377b8ee87 Mart*0231           PRESS0(i,j,bi,bj)=PRESS0(i,j,bi,bj)*HEFFM(i,j,bi,bj)
809c36b928 Patr*0232          ENDDO
                0233         ENDDO
                0234        ENDDO
                0235       ENDDO
                0236
                0237       IF ( SEAICEuseDYNAMICS ) THEN
                0238
7109a141b2 Patr*0239 #ifdef ALLOW_AUTODIFF_TAMC
                0240 CADJ STORE uice = comlev1, key=ikey_dynamics
                0241 CADJ STORE vice = comlev1, key=ikey_dynamics
                0242 #endif /* ALLOW_AUTODIFF_TAMC */
                0243
3f31c7d5de Mart*0244 crg what about ETA,ZETA
baa476eeba Dimi*0245
                0246 crg later c$taf loop = iteration uice,vice
                0247
09510da3bb Dimi*0248 cdm c$taf store uice,vice = comlev1_seaice_ds,
600feedf01 Dimi*0249 cdm c$taf&                key = kii + (ikey_dynamics-1)
809c36b928 Patr*0250 C NOW DO PREDICTOR TIME STEP
                0251       DO bj=myByLo(myThid),myByHi(myThid)
                0252        DO bi=myBxLo(myThid),myBxHi(myThid)
                0253         DO j=1-OLy,sNy+OLy
                0254          DO i=1-OLx,sNx+OLx
8377b8ee87 Mart*0255           UICENM1(i,j,bi,bj)=UICE(i,j,bi,bj)
                0256           VICENM1(i,j,bi,bj)=VICE(i,j,bi,bj)
45315406aa Mart*0257           UICEB(i,j,bi,bj)=UICE(i,j,bi,bj)
                0258           VICEB(i,j,bi,bj)=VICE(i,j,bi,bj)
809c36b928 Patr*0259          ENDDO
                0260         ENDDO
                0261        ENDDO
                0262       ENDDO
                0263
09510da3bb Dimi*0264 C NOW LSR SCHEME (ZHANG-J/HIBLER 1997)
7109a141b2 Patr*0265 CADJ STORE uice = comlev1, key=ikey_dynamics
                0266 CADJ STORE vice = comlev1, key=ikey_dynamics
                0267       CALL LSR( 1, myThid )
                0268 CADJ STORE uice = comlev1, key=ikey_dynamics
                0269 CADJ STORE vice = comlev1, key=ikey_dynamics
09510da3bb Dimi*0270
                0271 C NOW DO MODIFIED EULER STEP
                0272       DO bj=myByLo(myThid),myByHi(myThid)
                0273        DO bi=myBxLo(myThid),myBxHi(myThid)
                0274         DO j=1-OLy,sNy+OLy
                0275          DO i=1-OLx,sNx+OLx
8377b8ee87 Mart*0276           UICE(i,j,bi,bj)=HALF*(UICE(i,j,bi,bj)+UICENM1(i,j,bi,bj))
                0277           VICE(i,j,bi,bj)=HALF*(VICE(i,j,bi,bj)+VICENM1(i,j,bi,bj))
45315406aa Mart*0278           UICEB(i,j,bi,bj)=UICE(i,j,bi,bj)
                0279           VICEB(i,j,bi,bj)=VICE(i,j,bi,bj)
809c36b928 Patr*0280          ENDDO
                0281         ENDDO
                0282        ENDDO
                0283       ENDDO
                0284
09510da3bb Dimi*0285 C NOW LSR SCHEME (ZHANG-J/HIBLER 1997)
7109a141b2 Patr*0286       CALL LSR( 2, myThid )
809c36b928 Patr*0287
3f5db9535e Dimi*0288 cdm c$taf store uice,vice = comlev1, key=ikey_dynamics
baa476eeba Dimi*0289
809c36b928 Patr*0290       ENDIF
                0291
f9fa432bd6 Mart*0292 #ifdef ALLOW_DIAGNOSTICS
                0293       IF ( useDiagnostics ) THEN
                0294        CALL DIAGNOSTICS_FILL(zeta   ,'SIzeta  ',0,1 ,0,1,1,myThid)
                0295        CALL DIAGNOSTICS_FILL(eta    ,'SIeta   ',0,1 ,0,1,1,myThid)
                0296        CALL DIAGNOSTICS_FILL(press  ,'SIpress ',0,1 ,0,1,1,myThid)
                0297       ENDIF
5d6433c292 Jean*0298 #endif /* ALLOW_DIAGNOSTICS */
f9fa432bd6 Mart*0299
809c36b928 Patr*0300 C Calculate ocean surface stress
53092bcb42 Mart*0301       CALL OSTRES ( COR_ICE, myThid )
809c36b928 Patr*0302
d37c75ccbc Mart*0303 #ifdef SEAICE_ALLOW_CLIPVELS
                0304       IF ( SEAICEuseDYNAMICS .AND. SEAICE_clipVelocities) THEN
7109a141b2 Patr*0305 #ifdef ALLOW_AUTODIFF_TAMC
                0306 CADJ STORE uice = comlev1, key=ikey_dynamics
                0307 CADJ STORE vice = comlev1, key=ikey_dynamics
                0308 #endif /* ALLOW_AUTODIFF_TAMC */
809c36b928 Patr*0309 c Put a cap on ice velocity
                0310 c limit velocity to 0.40 m s-1 to avoid potential CFL violations
                0311 c in open water areas (drift of zero thickness ice)
                0312       DO bj=myByLo(myThid),myByHi(myThid)
                0313        DO bi=myBxLo(myThid),myBxHi(myThid)
                0314         DO j=1-OLy,sNy+OLy
                0315          DO i=1-OLx,sNx+OLx
                0316 #ifdef SEAICE_DEBUG
                0317 c          write(*,'(2i4,2i2,f7.1,7f12.3)')
                0318 c     &      i,j,bi,bj,UVM(I,J,bi,bj),amass(i,j,bi,bj)
                0319 c     &     ,gwatx(I,J,bi,bj),gwaty(i,j,bi,bj)
                0320 c     &     ,forcex(I,J,bi,bj),forcey(i,j,bi,bj)
772590b63c Mart*0321 c     &     ,uice(i,j,bi,bj)
                0322 c     &     ,vice(i,j,bi,bj)
6060ec2938 Dimi*0323 #endif /* SEAICE_DEBUG */
772590b63c Mart*0324           UICE(i,j,bi,bj)=min(UICE(i,j,bi,bj),0.40 _d +00)
                0325           VICE(i,j,bi,bj)=min(VICE(i,j,bi,bj),0.40 _d +00)
7109a141b2 Patr*0326          ENDDO
                0327         ENDDO
                0328        ENDDO
                0329       ENDDO
                0330 #ifdef ALLOW_AUTODIFF_TAMC
                0331 CADJ STORE uice = comlev1, key=ikey_dynamics
                0332 CADJ STORE vice = comlev1, key=ikey_dynamics
                0333 #endif /* ALLOW_AUTODIFF_TAMC */
                0334       DO bj=myByLo(myThid),myByHi(myThid)
                0335        DO bi=myBxLo(myThid),myBxHi(myThid)
                0336         DO j=1-OLy,sNy+OLy
                0337          DO i=1-OLx,sNx+OLx
772590b63c Mart*0338           UICE(i,j,bi,bj)=max(UICE(i,j,bi,bj),-0.40 _d +00)
                0339           VICE(i,j,bi,bj)=max(VICE(i,j,bi,bj),-0.40 _d +00)
809c36b928 Patr*0340          ENDDO
                0341         ENDDO
                0342        ENDDO
                0343       ENDDO
                0344
                0345       ENDIF
d37c75ccbc Mart*0346 #endif /* SEAICE_ALLOW_CLIPVELS */
45315406aa Mart*0347 #endif /* SEAICE_BGRID_DYNAMICS */
809c36b928 Patr*0348
                0349       RETURN
                0350       END