C $Header: /u/gcmpack/MITgcm/pkg/streamice/streamice_driving_stress_ppm.F,v 1.1 2013/06/12 21:30:22 dgoldberg Exp $ C $Name: $ #include "STREAMICE_OPTIONS.h" C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP SUBROUTINE STREAMICE_DRIVING_STRESS_PPM( 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 ! we are in an "active" cell IF (Gi.eq.1.AND..NOT.STREAMICE_EW_periodic) THEN ! western boundary - only one sided possible IF (STREAMICE_hmask(i+1,j,bi,bj).eq.1.0) THEN ! cell to east is active sx = (surf_el_streamice(i+1,j,bi,bj)- & surf_el_streamice(i,j,bi,bj))/dxC(i+1,j,bi,bj) ELSE ! cell to east is empty sx = 0. _d 0 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 * & (streamice_bg_surf_slope_x+sx) * & H_streamice(i,j,bi,bj) * rA(i,j,bi,bj) ENDIF ENDDO ENDDO ELSEIF (Gi.eq.Nx.AND..NOT.STREAMICE_EW_periodic) THEN ! eastern boundary - only one sided possible IF (STREAMICE_hmask(i-1,j,bi,bj).eq.1.0) THEN ! cell to west is active sx = (surf_el_streamice(i,j,bi,bj)- & surf_el_streamice(i-1,j,bi,bj))/dxC(i,j,bi,bj) ELSE ! cell to west is inactive sx = 0. _d 0 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 * & (streamice_bg_surf_slope_x+sx) * & H_streamice(i,j,bi,bj) * rA(i,j,bi,bj) ENDIF ENDDO ENDDO ELSE ! interior (west-east) cell IF (STREAMICE_hmask(i+1,j,bi,bj).eq.1.0 .and. & STREAMICE_hmask(i-1,j,bi,bj).eq.1.0) THEN k = 0 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) - & streamice_density * gravity * (1./6.) * & (-2.*surf_el_streamice(i-1,j,bi,bj) + & surf_el_streamice(i,j,bi,bj) + & surf_el_streamice(i+1,j,bi,bj) + & 3.*streamice_bg_surf_slope_x * dxF(i,j,bi,bj)) * & H_streamice(i,j,bi,bj) * .5 * dyF(i,j,bi,bj) ENDIF ENDDO k = 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) - & streamice_density * gravity * (1./6.) * & (-surf_el_streamice(i-1,j,bi,bj) - & surf_el_streamice(i,j,bi,bj) + & 2*surf_el_streamice(i+1,j,bi,bj) + & 3.*streamice_bg_surf_slope_x * dxF(i,j,bi,bj)) * & H_streamice(i,j,bi,bj) * .5 * dyF(i,j,bi,bj) ENDIF ENDDO ELSE 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) ELSEIF (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 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 * & (streamice_bg_surf_slope_x+sx) * & H_streamice(i,j,bi,bj) * rA(i,j,bi,bj) ENDIF ENDDO ENDDO ENDIF ENDIF !!!!!!!! DONE WITH X-GRADIENT IF (Gj.eq.1.AND..NOT.STREAMICE_NS_periodic) THEN ! western boundary - only one sided possible IF (STREAMICE_hmask(i,j+1,bi,bj).eq.1.0) THEN ! cell to east is active sy = (surf_el_streamice(i,j+1,bi,bj)- & surf_el_streamice(i,j,bi,bj))/dyC(i,j+1,bi,bj) ELSE ! cell to east is empty sy = 0. _d 0 ENDIF DO k=0,1 DO l=0,1 IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN taudy_SI(i+k,j+l,bi,bj) = taudy_SI(i+k,j+l,bi,bj) - & 0.25 * streamice_density * gravity * & (streamice_bg_surf_slope_y+sy) * & H_streamice(i,j,bi,bj) * rA(i,j,bi,bj) ENDIF ENDDO ENDDO ELSEIF (Gj.eq.Ny.AND..NOT.STREAMICE_NS_periodic) THEN ! eastern boundary - only one sided possible IF (STREAMICE_hmask(i,j-1,bi,bj).eq.1.0) THEN ! cell to west is active sy = (surf_el_streamice(i,j,bi,bj)- & surf_el_streamice(i,j-1,bi,bj))/dyC(i,j,bi,bj) ELSE ! cell to west is inactive sy = 0. _d 0 ENDIF DO k=0,1 DO l=0,1 IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN taudy_SI(i+k,j+l,bi,bj) = taudy_SI(i+k,j+l,bi,bj) - & 0.25 * streamice_density * gravity * & (streamice_bg_surf_slope_y+sy) * & H_streamice(i,j,bi,bj) * rA(i,j,bi,bj) ENDIF ENDDO ENDDO ELSE ! interior (west-east) cell IF (STREAMICE_hmask(i,j+1,bi,bj).eq.1.0 .and. & STREAMICE_hmask(i,j-1,bi,bj).eq.1.0) THEN l = 0 DO k=0,1 IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN taudy_SI(i+k,j+l,bi,bj) = taudy_SI(i+k,j+l,bi,bj) - & streamice_density * gravity * (1./6.) * & (-2.*surf_el_streamice(i,j-1,bi,bj) + & surf_el_streamice(i,j,bi,bj) + & surf_el_streamice(i,j+1,bi,bj) + & 3.*streamice_bg_surf_slope_y * dyF(i,j,bi,bj)) * & H_streamice(i,j,bi,bj) * .5 * dxF(i,j,bi,bj) ENDIF ENDDO l = 1 DO k=0,1 IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN taudy_SI(i+k,j+l,bi,bj) = taudy_SI(i+k,j+l,bi,bj) - & streamice_density * gravity * (1./6.) * & (-surf_el_streamice(i,j-1,bi,bj) - & surf_el_streamice(i,j,bi,bj) + & 2*surf_el_streamice(i,j+1,bi,bj) + & 3.*streamice_bg_surf_slope_y * dyF(i,j,bi,bj)) * & H_streamice(i,j,bi,bj) * .5 * dxF(i,j,bi,bj) ENDIF ENDDO ELSE 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))/dxC(i,j+1,bi,bj) ELSEIF (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))/dxC(i,j,bi,bj) ELSE sy = 0. _d 0 ENDIF DO k=0,1 DO l=0,1 IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN taudy_SI(i+k,j+l,bi,bj) = taudy_SI(i+k,j+l,bi,bj) - & 0.25 * streamice_density * gravity * & (streamice_bg_surf_slope_y+sy) * & H_streamice(i,j,bi,bj) * rA(i,j,bi,bj) ENDIF ENDDO ENDDO 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 * ! & (streamice_bg_surf_slope_x+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 * ! & (streamice_bg_surf_slope_y+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 #ifdef USE_ALT_RLOW neu_val = .5 * gravity * & (streamice_density * H_streamice (i,j,bi,bj) ** 2 - & streamice_density_ocean_avg * R_low_si(i,j,bi,bj) ** 2) #else neu_val = .5 * gravity * & (streamice_density * H_streamice (i,j,bi,bj) ** 2 - & streamice_density_ocean_avg * R_low(i,j,bi,bj) ** 2) #endif 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+streamice_addl_backstress) taudx_SI(i,j+1,bi,bj) = taudx_SI(i,j+1,bi,bj) - & .5 * dyG(i,j,bi,bj)*(neu_val+streamice_addl_backstress) 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+streamice_addl_backstress) ! 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+streamice_addl_backstress) 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+streamice_addl_backstress) taudy_SI(i+1,j,bi,bj) = taudy_SI(i+1,j,bi,bj) - & .5 * dxG(i,j,bi,bj)*(neu_val+streamice_addl_backstress) 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+streamice_addl_backstress) 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+streamice_addl_backstress) ENDIF ENDIF ENDDO ENDDO ENDDO ENDDO #endif RETURN END