pkg/seaice/seaice_bottomdrag_coeffs.F

9f4d410e86 Mart*0001 #include "SEAICE_OPTIONS.h"
                0002 #ifdef ALLOW_OBCS
                0003 # include "OBCS_OPTIONS.h"
                0004 #endif
                0005
                0006 CBOP
                0007 C     !ROUTINE: SEAICE_BOTTOMDRAG_COEFFS
                0008 C     !INTERFACE:
                0009       SUBROUTINE SEAICE_BOTTOMDRAG_COEFFS(
ec0d7df165 Mart*0010      I     uIceLoc, vIceLoc, HEFFMLoc,
9f4d410e86 Mart*0011 #ifdef SEAICE_ITD
ec0d7df165 Mart*0012      I     HEFFITDLoc, AREAITDLoc, AREALoc,
9f4d410e86 Mart*0013 #else
ec0d7df165 Mart*0014      I     HEFFLoc, AREALoc,
                0015 #endif
                0016      O     CbotLoc,
9f4d410e86 Mart*0017      I     iStep, myTime, myIter, myThid )
                0018
                0019 C     !DESCRIPTION: \bv
                0020 C     *==========================================================*
                0021 C     | SUBROUTINE SEAICE_BOTTOMDRAG_COEFFS
ec0d7df165 Mart*0022 C     | o Compute the non-linear drag coefficients for ice-bottom
9f4d410e86 Mart*0023 C     |   drag, as a parameterization for grounding fastice
ec0d7df165 Mart*0024 C     |   following
9f4d410e86 Mart*0025 C     |   Lemieux et al. (2015), doi:10.1002/2014JC010678
                0026 C     *==========================================================*
                0027 C     | written by Martin Losch, Apr 2016
                0028 C     *==========================================================*
                0029 C     \ev
                0030
                0031 C     !USES:
                0032       IMPLICIT NONE
                0033
                0034 C     === Global variables ===
                0035 #include "SIZE.h"
                0036 #include "EEPARAMS.h"
                0037 #include "PARAMS.h"
                0038 #include "GRID.h"
                0039 #include "DYNVARS.h"
                0040 #include "SEAICE_SIZE.h"
                0041 #include "SEAICE_PARAMS.h"
                0042
                0043 C     !INPUT/OUTPUT PARAMETERS:
                0044 C     === Routine arguments ===
ec0d7df165 Mart*0045 C     u/vIceLoc :: local copies of the current ice velocity
                0046 C     HEFFMLoc  :: local copy of land-sea masks
                0047 C     CbotLoc   :: drag coefficients
                0048 C     iStep     :: current sub-time step iterate
                0049 C     myTime    :: Simulation time
                0050 C     myIter    :: Simulation timestep number
                0051 C     myThid    :: my Thread Id. number
d5254b4e3d Mart*0052       _RL uIceLoc   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0053       _RL vIceLoc   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
ec0d7df165 Mart*0054       _RL HEFFMLoc  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0055 #ifdef SEAICE_ITD
                0056       _RL HEFFITDLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nITD,nSx,nSy)
                0057       _RL AREAITDLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nITD,nSx,nSy)
                0058 #else
d5254b4e3d Mart*0059       _RL HEFFLoc   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
ec0d7df165 Mart*0060 #endif
d5254b4e3d Mart*0061       _RL AREALoc   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0062       _RL CbotLoc   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
ec0d7df165 Mart*0063       INTEGER iStep
9f4d410e86 Mart*0064       _RL     myTime
                0065       INTEGER myIter
                0066       INTEGER myThid
                0067
                0068 #ifdef SEAICE_ALLOW_BOTTOMDRAG
                0069 C     === local variables ===
                0070 C     i,j,bi,bj,ksrf :: loop indices
                0071       INTEGER i,j,bi,bj
                0072       INTEGER kSrf
                0073 #ifdef SEAICE_ITD
                0074       INTEGER k
                0075 #endif /* SEAICE_ITD */
                0076       _RL     tmpFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0077       _RL     tmp, hActual, hCrit, recip_k1, u0sq, fac, rFac
                0078 CEOP
                0079
d5254b4e3d Mart*0080       IF (SEAICEbasalDragK2.GT.0. _d 0) THEN
                0081 C     avoid this computation for a non-zero coefficient
0320e25227 Mart*0082       IF ( usingPCoords ) THEN
                0083        kSrf = Nr
                0084       ELSE
                0085        kSrf = 1
                0086       ENDIF
9f4d410e86 Mart*0087 C     some abbreviations
                0088       u0sq     = SEAICEbasalDragU0*SEAICEbasalDragU0
                0089       recip_k1 = 0. _d 0
ec0d7df165 Mart*0090       IF ( SEAICEbasalDragK1 .GT. 0. _d 0 )
9f4d410e86 Mart*0091      &     recip_k1 = 1. _d 0/SEAICEbasalDragK1
                0092 C     fac scales the soft maximum for more accuracy
                0093       fac = 10. _d 0
                0094       rFac = 1./fac
                0095
                0096       DO bj=myByLo(myThid),myByHi(myThid)
                0097        DO bi=myBxLo(myThid),myBxHi(myThid)
ec0d7df165 Mart*0098         DO j=1-OLy,sNy+OLy
4a08d54d3a Mart*0099          DO i=1-OLx,sNx+OLx
d5254b4e3d Mart*0100           CbotLoc(i,j,bi,bj) = 0. _d 0
                0101           tmpFld (i,j)       = 0. _d 0
9f4d410e86 Mart*0102          ENDDO
                0103         ENDDO
ec0d7df165 Mart*0104         DO j=1-OLy,sNy+OLy-1
4a08d54d3a Mart*0105          DO i=1-OLx,sNx+OLx-1
d5254b4e3d Mart*0106           IF ( AREALoc(i,j,bi,bj) .GT. 0.01 _d 0 ) THEN
9f4d410e86 Mart*0107 #ifdef OBCS_UVICE_OLD
                0108            tmp = 0.25 _d 0*(
ec0d7df165 Mart*0109      &          (   uIceLoc(i  ,j,bi,bj)+uIceLoc(i+1,j,bi,bj)
9f4d410e86 Mart*0110      &          )**2
ec0d7df165 Mart*0111      &          + ( vIceLoc(i, j ,bi,bj)+vIceLoc(i,j+1,bi,bj)
9f4d410e86 Mart*0112      &          )**2 )
                0113 #else /* OBCS_UVICE_OLD */
                0114            tmp = 0.25 _d 0*(
ec0d7df165 Mart*0115      &          ( uIceLoc(i  ,j,bi,bj)*maskInW( i ,j,bi,bj)
d5254b4e3d Mart*0116      &          + uIceLoc(i+1,j,bi,bj)*maskInW(i+1,j,bi,bj) )**2
                0117      &        + ( vIceLoc(i,j  ,bi,bj)*maskInS(i,j  ,bi,bj)
ec0d7df165 Mart*0118      &          + vIceLoc(i,j+1,bi,bj)*maskInS(i,j+1,bi,bj) )**2 )
9f4d410e86 Mart*0119 #endif /* OBCS_UVICE_OLD */
d5254b4e3d Mart*0120            tmpFld(i,j) = SEAICEbasalDragK2 / SQRT(tmp + u0sq)
9f4d410e86 Mart*0121           ENDIF
                0122          ENDDO
                0123         ENDDO
                0124 #ifdef SEAICE_ITD
                0125         DO k=1,nITD
                0126 #endif /* SEAICE_ITD */
ec0d7df165 Mart*0127          DO j=1-OLy,sNy+OLy-1
4a08d54d3a Mart*0128           DO i=1-OLx,sNx+OLx-1
d5254b4e3d Mart*0129            IF ( AREALoc(i,j,bi,bj) .GT. 0.01 _d 0 ) THEN
ec0d7df165 Mart*0130 CML           hActual = HEFFLoc(i,j,bi,bj)
                0131 CML     &          /SQRT( AREAITDLoc(i,j,bi,bj)**2 + area_reg_sq )
9f4d410e86 Mart*0132 CML           hActual = SQRT(hActual * hActual + hice_reg_sq)
d5254b4e3d Mart*0133 CML           hCrit   = ABS(R_low(i,j,bi,bj)) * recip_k1
9f4d410e86 Mart*0134 #ifdef SEAICE_ITD
d5254b4e3d Mart*0135             hActual = HEFFITDLoc(i,j,k,bi,bj)
9f4d410e86 Mart*0136 C     here we do not need recip_k1, because we resolve the very thick ice
d5254b4e3d Mart*0137             hCrit   = ABS(R_low(i,j,bi,bj))*AREAITDLoc(i,j,k,bi,bj)
9f4d410e86 Mart*0138 #else
d5254b4e3d Mart*0139             hActual = HEFFLoc(i,j,bi,bj)
                0140             hCrit   = ABS(R_low(i,j,bi,bj))*AREALoc(i,j,bi,bj)*recip_k1
9f4d410e86 Mart*0141 #endif /* SEAICE_ITD */
                0142 C     we want to have some soft maximum for better differentiability:
                0143 C     max(a,b;k) = ln(exp(k*a)+exp(k*b))/k
                0144 C     In our case, b=0, so exp(k*b) = 1.
                0145 C     max(a,0;k) = ln(exp(k*a)+1)/k
                0146 C     If k*a gets too large, EXP will overflow, but for the anticipated
                0147 C     values of hActual < 100m, and k=10, this should be very unlikely
d5254b4e3d Mart*0148 CML             CbotLoc(i,j,bi,bj) =
                0149 CML     &            tmpFld(i,j) * MAX( hActual - hCrit, 0. _d 0)
                0150             CbotLoc(i,j,bi,bj) = CbotLoc(i,j,bi,bj)
                0151      &           + tmpFld(i,j)
9f4d410e86 Mart*0152      &           * LOG(EXP( fac*(hActual - hCrit) ) + 1. _d 0)*rFac
                0153      &           * EXP( - SEAICE_cBasalStar
d5254b4e3d Mart*0154      &                  *(SEAICE_area_max - AREALoc(i,j,bi,bj)) )
                0155      &           * HEFFMLoc(i,j,bi,bj)
9f4d410e86 Mart*0156            ENDIF
                0157           ENDDO
                0158          ENDDO
                0159 #ifdef SEAICE_ITD
                0160         ENDDO
                0161 #endif /* SEAICE_ITD */
                0162        ENDDO
                0163       ENDDO
d5254b4e3d Mart*0164 C     endif SEAICEbasalDragK2.GT.0.
                0165       ENDIF
9f4d410e86 Mart*0166 #endif /* SEAICE_ALLOW_BOTTOMDRAG */
                0167
                0168       RETURN
                0169       END