C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_bottomdrag_coeffs.F,v 1.1 2016/04/22 08:43:41 mlosch Exp $ C $Name: $ #include "SEAICE_OPTIONS.h" #ifdef ALLOW_OBCS # include "OBCS_OPTIONS.h" #endif #ifdef ALLOW_AUTODIFF # include "AUTODIFF_OPTIONS.h" #endif CBOP C !ROUTINE: SEAICE_BOTTOMDRAG_COEFFS C !INTERFACE: SUBROUTINE SEAICE_BOTTOMDRAG_COEFFS( I uIce, vIce, #ifdef SEAICE_ITD I HEFFITD, AREAITD, AREA, #else I HEFF, AREA, #endif O CbotC, I iStep, myTime, myIter, myThid ) C !DESCRIPTION: \bv C *==========================================================* C | SUBROUTINE SEAICE_BOTTOMDRAG_COEFFS C | o Compute the non-linear drag coefficients for ice-bottom C | drag, as a parameterization for grounding fastice C | following C | Lemieux et al. (2015), doi:10.1002/2014JC010678 C *==========================================================* C | written by Martin Losch, Apr 2016 C *==========================================================* C \ev C !USES: IMPLICIT NONE C === Global variables === #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" #include "DYNVARS.h" #include "SEAICE_SIZE.h" #include "SEAICE_PARAMS.h" C !INPUT/OUTPUT PARAMETERS: C === Routine arguments === C myTime :: Simulation time C myIter :: Simulation timestep number C myThid :: my Thread Id. number C iStep :: current sub-time step iterate _RL myTime INTEGER myIter INTEGER myThid INTEGER iStep C u/vIce :: local copies of the current ice velocity _RL uIce(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL vIce(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) #ifdef SEAICE_ITD _RL HEFFITD(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nITD,nSx,nSy) _RL AREAITD(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nITD,nSx,nSy) #else _RL HEFF(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) #endif _RL AREA(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C CbotC :: drag coefficients _RL CbotC(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) #ifdef SEAICE_ALLOW_BOTTOMDRAG C === local variables === C i,j,bi,bj,ksrf :: loop indices INTEGER i,j,bi,bj INTEGER kSrf #ifdef SEAICE_ITD INTEGER k #endif /* SEAICE_ITD */ _RL tmpFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL tmp, hActual, hCrit, recip_k1, u0sq, fac, rFac CEOP kSrf=1 C some abbreviations u0sq = SEAICEbasalDragU0*SEAICEbasalDragU0 recip_k1 = 0. _d 0 IF ( SEAICEbasalDragK1 .GT. 0. _d 0 ) & recip_k1 = 1. _d 0/SEAICEbasalDragK1 C fac scales the soft maximum for more accuracy fac = 10. _d 0 rFac = 1./fac DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1-OLy,sNy+Oly DO i=1-OLx,sNx+OLy CbotC(I,J,bi,bj) = 0. _d 0 tmpFld(I,J) = 0. _d 0 ENDDO ENDDO DO j=1-OLy,sNy+Oly-1 DO i=1-OLx,sNx+OLy-1 IF ( AREA(I,J,bi,bj) .GT. 0.01 _d 0 ) THEN #ifdef OBCS_UVICE_OLD tmp = 0.25 _d 0*( & ( uIce(i ,j,bi,bj)+uIce(i+1,j,bi,bj) & )**2 & + ( vIce(i, j ,bi,bj)+vIce(i,j+1,bi,bj) & )**2 ) #else /* OBCS_UVICE_OLD */ tmp = 0.25 _d 0*( & ( uIce(i ,j,bi,bj)*maskInW( i ,j,bi,bj) & + uIce(i+1,J,bi,bj)*maskInW(i+1,j,bi,bj) )**2 & + ( vIce(i,j ,bi,bj)*maskInS(i, j ,bi,bj) & + vIce(i,j+1,bi,bj)*maskInS(i,j+1,bi,bj) )**2 ) #endif /* OBCS_UVICE_OLD */ tmpFld(I,J) = SEAICEbasalDragK2 / SQRT(tmp + u0sq) ENDIF ENDDO ENDDO #ifdef SEAICE_ITD DO k=1,nITD #endif /* SEAICE_ITD */ DO j=1-OLy,sNy+Oly-1 DO i=1-OLx,sNx+OLy-1 IF ( AREA(I,J,bi,bj) .GT. 0.01 _d 0 ) THEN CML hActual = HEFF(i,j,bi,bj) CML & /SQRT( AREAITD(i,j,bi,bj)**2 + area_reg_sq ) CML hActual = SQRT(hActual * hActual + hice_reg_sq) CML hCrit = ABS(R_low(I,J,bi,bj)) * recip_k1 #ifdef SEAICE_ITD hActual = HEFFITD(I,J,k,bi,bj) C here we do not need recip_k1, because we resolve the very thick ice hCrit = ABS(R_low(I,J,bi,bj))*AREAITD(I,J,k,bi,bj) #else hActual = HEFF(I,J,bi,bj) hCrit = ABS(R_low(I,J,bi,bj))*AREA(I,J,bi,bj)*recip_k1 #endif /* SEAICE_ITD */ C we want to have some soft maximum for better differentiability: C max(a,b;k) = ln(exp(k*a)+exp(k*b))/k C In our case, b=0, so exp(k*b) = 1. C max(a,0;k) = ln(exp(k*a)+1)/k C If k*a gets too large, EXP will overflow, but for the anticipated C values of hActual < 100m, and k=10, this should be very unlikely CML CbotC(I,J,bi,bj) = CML & tmpFld(I,J) * MAX( hActual - hCrit, 0. _d 0) CbotC(I,J,bi,bj) = CbotC(I,J,bi,bj) & + tmpFld(I,J) & * LOG(EXP( fac*(hActual - hCrit) ) + 1. _d 0)*rFac & * EXP( - SEAICE_cBasalStar & *(SEAICE_area_max - AREA(I,J,bi,bj)) ) & * maskC(I,J,kSrf,bi,bj) ENDIF ENDDO ENDDO #ifdef SEAICE_ITD ENDDO #endif /* SEAICE_ITD */ ENDDO ENDDO #endif /* SEAICE_ALLOW_BOTTOMDRAG */ RETURN END