pkg/seaice/seaice_freedrift.F

e97d1f6d02 Gael*0001 #include "SEAICE_OPTIONS.h"
                0002
                0003 CBOP
                0004       SUBROUTINE SEAICE_FREEDRIFT( myTime, myIter, myThid )
850adf2a7e Jean*0005 C     *==========================================================*
                0006 C     | SUBROUTINE  SEAICE_FREEDRIFT
                0007 C     | o Solve ice approximate momentum equation analytically
                0008 C     *==========================================================*
e97d1f6d02 Gael*0009       IMPLICIT NONE
                0010
                0011 C     === Global variables ===
                0012 #include "SIZE.h"
                0013 #include "EEPARAMS.h"
                0014 #include "PARAMS.h"
                0015 #include "DYNVARS.h"
                0016 #include "GRID.h"
03c669d1ab Jean*0017 #include "SEAICE_SIZE.h"
e97d1f6d02 Gael*0018 #include "SEAICE_PARAMS.h"
3f0f10fc37 Mart*0019 #include "SEAICE_GRID.h"
03c669d1ab Jean*0020 #include "SEAICE.h"
e97d1f6d02 Gael*0021
                0022 C     === Routine arguments ===
                0023 C     myTime :: Simulation time
                0024 C     myIter :: Simulation timestep number
                0025 C     myThid :: my Thread Id. number
                0026       _RL     myTime
                0027       INTEGER myIter
                0028       INTEGER myThid
                0029 CEOP
                0030
45315406aa Mart*0031 #if ( defined SEAICE_CGRID && defined SEAICE_ALLOW_FREEDRIFT )
e97d1f6d02 Gael*0032
                0033 C     === Local variables ===
73c8571595 Jean*0034       INTEGER i, j, kSrf, bi, bj
e97d1f6d02 Gael*0035
                0036       _RL tmpscal1,tmpscal2,tmpscal3,tmpscal4
                0037
                0038       _RL taux_onIce_cntr, tauy_onIce_cntr, uvel_cntr, vvel_cntr
                0039       _RL mIceCor, rhs_x, rhs_y, rhs_n, rhs_a, sol_n, sol_a
                0040
850adf2a7e Jean*0041       _RL uice_cntr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0042       _RL vice_cntr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
e97d1f6d02 Gael*0043
0320e25227 Mart*0044       IF ( usingPCoords ) THEN
                0045        kSrf = Nr
                0046       ELSE
                0047        kSrf = 1
                0048       ENDIF
e97d1f6d02 Gael*0049
850adf2a7e Jean*0050 C initialize fields:
                0051 C ==================
e97d1f6d02 Gael*0052
0320e25227 Mart*0053       DO bj=myByLo(myThid),myByHi(myThid)
                0054        DO bi=myBxLo(myThid),myBxHi(myThid)
                0055         DO j=1-OLy,sNy+OLy
                0056          DO i=1-OLx,sNx+OLx
                0057           uice_fd(i,j,bi,bj)=0. _d 0
                0058           vice_fd(i,j,bi,bj)=0. _d 0
                0059           uice_cntr(i,j,bi,bj)=0. _d 0
                0060           Vice_cntr(i,j,bi,bj)=0. _d 0
e97d1f6d02 Gael*0061          ENDDO
                0062         ENDDO
                0063        ENDDO
0320e25227 Mart*0064       ENDDO
e97d1f6d02 Gael*0065
0320e25227 Mart*0066       DO bj=myByLo(myThid),myByHi(myThid)
                0067        DO bi=myBxLo(myThid),myBxHi(myThid)
                0068         DO j=1,sNy
                0069          DO i=1,sNx
73c8571595 Jean*0070
850adf2a7e Jean*0071 C preliminary computations:
                0072 C =========================
                0073 C air-ice stress at cell center
0320e25227 Mart*0074           taux_onIce_cntr=HALF*
                0075      &      (FORCEX0(i,j,bi,bj)+FORCEX0(i+1,j,bi,bj))
                0076           tauy_onIce_cntr=HALF*
                0077      &      (FORCEY0(i,j,bi,bj)+FORCEY0(i,j+1,bi,bj))
850adf2a7e Jean*0078 C mass of ice per unit area (kg/m2) times coriolis f
0320e25227 Mart*0079           mIceCor=SEAICE_rhoIce*HEFF(i,j,bi,bj)*_fCori(i,j,bi,bj)
850adf2a7e Jean*0080 C ocean velocity at cell center
0320e25227 Mart*0081           uvel_cntr=HALF*(uvel(i,j,kSrf,bi,bj)+uvel(i+1,j,kSrf,bi,bj))
                0082           vvel_cntr=HALF*(vvel(i,j,kSrf,bi,bj)+vvel(i,j+1,kSrf,bi,bj))
850adf2a7e Jean*0083 C right hand side of free drift equation:
0320e25227 Mart*0084           rhs_x= -taux_onIce_cntr -mIceCor*vvel_cntr
                0085           rhs_y= -tauy_onIce_cntr +mIceCor*uvel_cntr
e97d1f6d02 Gael*0086
850adf2a7e Jean*0087 C norm of angle of rhs
0320e25227 Mart*0088           tmpscal1=rhs_x*rhs_x + rhs_y*rhs_y
                0089           IF ( tmpscal1.GT.ZERO ) THEN
                0090            rhs_n=SQRT( rhs_x*rhs_x + rhs_y*rhs_y )
                0091            rhs_a=ATAN2(rhs_y,rhs_x)
                0092           ELSE
                0093            rhs_n=0. _d 0
                0094            rhs_a=0. _d 0
                0095           ENDIF
e97d1f6d02 Gael*0096
850adf2a7e Jean*0097 C solve for norm:
                0098 C ===============
0320e25227 Mart*0099           IF ( YC(i,j,bi,bj) .LT. ZERO ) THEN
                0100            tmpscal1 = 1. _d 0 /rhoConst/SEAICE_waterDrag_south
                0101           ELSE
                0102            tmpscal1 = 1. _d 0 /rhoConst/SEAICE_waterDrag
                0103           ENDIF
850adf2a7e Jean*0104 C polynomial coefficients
0320e25227 Mart*0105           tmpscal2= tmpscal1*tmpscal1*mIceCor*mIceCor
                0106           tmpscal3= tmpscal1*tmpscal1*rhs_n*rhs_n
850adf2a7e Jean*0107 C discriminant
0320e25227 Mart*0108           tmpscal4=tmpscal2*tmpscal2+4. _d 0*tmpscal3
                0109           IF ( tmpscal3.GT.ZERO ) THEN
                0110            sol_n=SQRT(HALF*(SQRT(tmpscal4)-tmpscal2))
                0111           ELSE
                0112            sol_n=0. _d 0
                0113           ENDIF
e97d1f6d02 Gael*0114
850adf2a7e Jean*0115 C solve for angle:
                0116 C ================
0320e25227 Mart*0117           IF ( YC(i,j,bi,bj) .LT. ZERO ) THEN
                0118            tmpscal1 = SEAICE_waterDrag_south*rhoConst
                0119           ELSE
                0120            tmpscal1 = SEAICE_waterDrag*rhoConst
                0121           ENDIF
                0122
                0123           tmpscal2= tmpscal1*sol_n*sol_n
                0124           tmpscal3= mIceCor*sol_n
                0125
                0126           tmpscal4=tmpscal2*tmpscal2 + tmpscal3*tmpscal3
                0127           IF ( tmpscal4.GT.ZERO ) THEN
                0128            sol_a=rhs_a-ATAN2(tmpscal3,tmpscal2)
                0129           ELSE
                0130            sol_a=0. _d 0
                0131           ENDIF
e97d1f6d02 Gael*0132
850adf2a7e Jean*0133 C compute uice, vice at cell center:
                0134 C ==================================
0320e25227 Mart*0135           uice_cntr(i,j,bi,bj)=uvel_cntr-sol_n*COS(sol_a)
                0136           vice_cntr(i,j,bi,bj)=vvel_cntr-sol_n*SIN(sol_a)
e97d1f6d02 Gael*0137
                0138          ENDDO
                0139         ENDDO
                0140        ENDDO
0320e25227 Mart*0141       ENDDO
e97d1f6d02 Gael*0142
850adf2a7e Jean*0143 C interpolated to velocity points:
                0144 C ================================
e97d1f6d02 Gael*0145
0320e25227 Mart*0146       CALL EXCH_UV_AGRID_3D_RL(uice_cntr,vice_cntr,.TRUE.,1,myThid)
                0147
                0148       DO bj=myByLo(myThid),myByHi(myThid)
                0149        DO bi=myBxLo(myThid),myBxHi(myThid)
                0150         DO j=1,sNy
                0151          DO i=1,sNx
                0152           uice_fd(i,j,bi,bj)=HALF*
                0153      &      (uice_cntr(i-1,j,bi,bj)+uice_cntr(i,j,bi,bj))
                0154           vice_fd(i,j,bi,bj)=HALF*
                0155      &      (vice_cntr(i,j-1,bi,bj)+vice_cntr(i,j,bi,bj))
e97d1f6d02 Gael*0156          ENDDO
                0157         ENDDO
                0158        ENDDO
0320e25227 Mart*0159       ENDDO
e97d1f6d02 Gael*0160
0320e25227 Mart*0161       CALL EXCH_UV_XY_RL( uice_fd, vice_fd, .TRUE., myThid )
e97d1f6d02 Gael*0162
73c8571595 Jean*0163 C     Apply masks (same/similar to seaice_evp.F/seaice_lsr.F)
                0164       DO bj=myByLo(myThid),myByHi(myThid)
                0165        DO bi=myBxLo(myThid),myBxHi(myThid)
850adf2a7e Jean*0166         DO j=1-OLy,sNy+OLy
                0167          DO i=1-OLx,sNx+OLx
ec0d7df165 Mart*0168           uIce_fd(i,j,bi,bj)=uIce_fd(i,j,bi,bj)*SIMaskU(i,j,bi,bj)
                0169           vIce_fd(i,j,bi,bj)=vIce_fd(i,j,bi,bj)*SIMaskV(i,j,bi,bj)
73c8571595 Jean*0170          ENDDO
                0171         ENDDO
                0172        ENDDO
                0173       ENDDO
                0174
45315406aa Mart*0175 #endif /* SEAICE_CGRID and SEAICE_ALLOW_FREEDRIFT */
e97d1f6d02 Gael*0176       RETURN
                0177       END