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