C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/dgoldberg/streamice/streamice_driving_stress.F,v 1.1 2012/03/29 15:59:21 heimbach Exp $ C $Name: $ #include "STREAMICE_OPTIONS.h" C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP SUBROUTINE STREAMICE_DRIVING_STRESS( myThid ) ! O taudx, ! O taudy ) C /============================================================\ C | SUBROUTINE | C | o | C |============================================================| C | | C \============================================================/ IMPLICIT NONE C === Global variables === #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" #include "STREAMICE.h" #include "STREAMICE_CG.h" C !INPUT/OUTPUT ARGUMENTS INTEGER myThid ! _RL taudx (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) ! _RL taudx (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) #ifdef ALLOW_STREAMICE C LOCAL VARIABLES INTEGER i, j, bi, bj, k, l, & Gi, Gj LOGICAL at_west_bdry, at_east_bdry, & at_north_bdry, at_south_bdry _RL sx, sy, diffx, diffy, neu_val IF (myXGlobalLo.eq.1) at_west_bdry = .true. IF (myYGlobalLo.eq.1) at_south_bdry = .true. IF (myXGlobalLo-1+sNx*nSx.eq.Nx) & at_east_bdry = .false. IF (myYGlobalLo-1+sNy*nSy.eq.Ny) & at_north_bdry = .false. DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx taudx_SI(i,j,bi,bj) = 0. _d 0 taudy_SI(i,j,bi,bj) = 0. _d 0 ENDDO ENDDO ENDDO ENDDO DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO i=0,sNx+1 DO j=0,sNy+1 diffx = 0. _d 0 diffy = 0. _d 0 sx = 0. _d 0 sy = 0. _d 0 Gi = (myXGlobalLo-1)+(bi-1)*sNx+i Gj = (myYGlobalLo-1)+(bj-1)*sNy+j IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN IF (Gi .eq. 1) THEN IF (STREAMICE_hmask(i+1,j,bi,bj).eq.1.0) THEN sx = (surf_el_streamice(i+1,j,bi,bj)- & surf_el_streamice(i,j,bi,bj))/dxC(i+1,j,bi,bj) ELSE sx = 0. _d 0 ENDIF ELSEIF (Gi .eq. Nx) THEN IF (STREAMICE_hmask(i-1,j,bi,bj).eq.1.0) THEN sx = (surf_el_streamice(i,j,bi,bj)- & surf_el_streamice(i-1,j,bi,bj))/dxC(i,j,bi,bj) ELSE sx = 0. _d 0 ENDIF ELSE IF (STREAMICE_hmask(i+1,j,bi,bj).eq.1.0) THEN diffx = diffx + dxC(i+1,j,bi,bj) sx = surf_el_streamice(i+1,j,bi,bj) ELSE sx = surf_el_streamice(i,j,bi,bj) ENDIF IF (STREAMICE_hmask(i-1,j,bi,bj).eq.1.0) THEN diffx = diffx + dxC(i,j,bi,bj) sx = sx - surf_el_streamice(i-1,j,bi,bj) ELSE sx = sx - surf_el_streamice(i,j,bi,bj) ENDIF IF (diffx .gt. 0. _d 0) THEN sx = sx / diffx ELSE sx = 0. _d 0 ENDIF ENDIF IF (Gj .eq. 1) THEN IF (STREAMICE_hmask(i,j+1,bi,bj).eq.1.0) THEN sy = (surf_el_streamice(i,j+1,bi,bj)- & surf_el_streamice(i,j,bi,bj))/dyC(i,j+1,bi,bj) ELSE sy = 0. _d 0 ENDIF ELSEIF (Gj .eq. Ny) THEN IF (STREAMICE_hmask(i,j-1,bi,bj).eq.1.0) THEN sy = (surf_el_streamice(i,j,bi,bj)- & surf_el_streamice(i,j-1,bi,bj))/dyC(i,j,bi,bj) ELSE sy = 0. _d 0 ENDIF ELSE IF (STREAMICE_hmask(i,j+1,bi,bj).eq.1.0) THEN diffy = diffy + dyC(i,j+1,bi,bj) sy = surf_el_streamice(i,j+1,bi,bj) ELSE sy = surf_el_streamice(i,j,bi,bj) ENDIF IF (STREAMICE_hmask(i,j-1,bi,bj).eq.1.0) THEN diffy = diffy + dyC(i,j,bi,bj) sy = sy - surf_el_streamice(i,j-1,bi,bj) ELSE sy = sy - surf_el_streamice(i,j,bi,bj) ENDIF IF (diffy .gt. 0. _d 0) THEN sy = sy / diffy ELSE sy = 0. _d 0 ENDIF ENDIF DO k=0,1 DO l=0,1 IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN taudx_SI(i+k,j+l,bi,bj) = taudx_SI(i+k,j+l,bi,bj) - & 0.25 * streamice_density * gravity * sx * & H_streamice(i,j,bi,bj) * rA(i,j,bi,bj) taudy_SI(i+k,j+l,bi,bj) = taudy_SI(i+k,j+l,bi,bj) - & 0.25 * streamice_density * gravity * sy * & H_streamice(i,j,bi,bj) * rA(i,j,bi,bj) ENDIF ENDDO ENDDO IF (float_frac_streamice(i,j,bi,bj) .eq. 1.0) then neu_val = .5 * gravity * & (streamice_density * H_streamice (i,j,bi,bj) ** 2 - & streamice_density_ocean_avg * R_low(i,j,bi,bj) ** 2) ELSE neu_val = .5 * gravity * & (1-streamice_density/streamice_density_ocean_avg) * & streamice_density * H_streamice(i,j,bi,bj) ** 2 ENDIF IF ((STREAMICE_ufacemask(i,j,bi,bj) .eq. 2) & .OR. (STREAMICE_hmask(i-1,j,bi,bj) .eq. 0) & .OR. (STREAMICE_hmask(i-1,j,bi,bj) .eq. 2) ) THEN ! left face of the cell is at a stress boundary ! the depth-integrated longitudinal stress is equal to the difference of depth-integrated pressure on either side of the face ! on the ice side, it is rho g h^2 / 2 ! on the ocean side, it is rhow g (delta OD)^2 / 2 ! OD can be zero under the ice; but it is ASSUMED on the ice-free side of the face, topography elevation is not above the base of the ! ice in the current cell taudx_SI(i,j,bi,bj) = taudx_SI(i,j,bi,bj) - & .5 * dyG(i,j,bi,bj) * neu_val ! note negative sign is due to direction of normal vector taudx_SI(i,j+1,bi,bj) = taudx_SI(i,j+1,bi,bj) - & .5 * dyG(i,j,bi,bj) * neu_val ENDIF IF ((STREAMICE_ufacemask(i+1,j,bi,bj) .eq. 2) & .OR. (STREAMICE_hmask(i+1,j,bi,bj) .eq. 0) & .OR. (STREAMICE_hmask(i+1,j,bi,bj) .eq. 2) ) THEN taudx_SI(i+1,j,bi,bj) = taudx_SI(i+1,j,bi,bj) + & .5 * dyG(i+1,j,bi,bj) * neu_val ! note negative sign is due to direction of normal vector taudx_SI(i+1,j+1,bi,bj) = taudx_SI(i+1,j+1,bi,bj) + & .5 * dyG(i+1,j,bi,bj) * neu_val ENDIF IF ((STREAMICE_vfacemask(i,j,bi,bj) .eq. 2) & .OR. (STREAMICE_hmask(i,j-1,bi,bj) .eq. 0) & .OR. (STREAMICE_hmask(i,j-1,bi,bj) .eq. 2) ) THEN taudy_SI(i,j,bi,bj) = taudy_SI(i,j,bi,bj) - & .5 * dxG(i,j,bi,bj) * neu_val ! note negative sign is due to direction of normal vector taudy_SI(i+1,j,bi,bj) = taudy_SI(i+1,j,bi,bj) - & .5 * dxG(i,j,bi,bj) * neu_val ENDIF IF ((STREAMICE_vfacemask(i,j+1,bi,bj) .eq. 2) & .OR. (STREAMICE_hmask(i,j+1,bi,bj) .eq. 0) & .OR. (STREAMICE_hmask(i,j+1,bi,bj) .eq. 2) ) THEN taudy_SI(i,j+1,bi,bj) = taudy_SI(i,j+1,bi,bj) + & .5 * dxG(i,j+1,bi,bj) * neu_val ! note negative sign is due to direction of normal vector taudy_SI(i+1,j+1,bi,bj) = taudy_SI(i+1,j+1,bi,bj) + & .5 * dxG(i,j+1,bi,bj) * neu_val ENDIF ENDIF ENDDO ENDDO ENDDO ENDDO #endif RETURN END