C $Header: /u/gcmpack/MITgcm/pkg/streamice/streamice_driving_stress.F,v 1.10 2015/10/18 14:37:49 dgoldberg 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 grd_below_sl
_RL sx, sy, diffx, diffy, geom_fac
_RL i_rhow
_RL avg_density
& (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
#ifdef STREAMICE_FIRN_CORRECTION
_RL firn_depth, h
#endif
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.
#ifdef STREAMICE_FIRN_CORRECTION
firn_depth = streamice_density *
& streamice_firn_correction
& / (streamice_density-streamice_density_firn)
#endif
i_rhow = 1./streamice_density_ocean_avg
DO bj = myByLo(myThid), myByHi(myThid)
DO bi = myBxLo(myThid), myBxHi(myThid)
DO j=1-OLy+1,sNy+OLy-1
DO i=1-OLx+1,sNx+OLx-1
taudx_SI(i,j,bi,bj) = 0. _d 0
taudy_SI(i,j,bi,bj) = 0. _d 0
#ifdef STREAMICE_FIRN_CORRECTION
if (STREAMICE_apply_firn_correction) then
if (streamice_hmask(i,j,bi,bj).eq.1) then
h = h_streamice(i,j,bi,bj)
if (h.lt.firn_depth) then
avg_density(i,j,bi,bj) = streamice_density_firn
else
avg_density(i,j,bi,bj) = streamice_density *
& (h - streamice_firn_correction) / h
endif
endif
else
#endif
avg_density(i,j,bi,bj) = streamice_density
#ifdef STREAMICE_FIRN_CORRECTION
endif
#endif
ENDDO
ENDDO
ENDDO
ENDDO
DO bj = myByLo(myThid), myByHi(myThid)
DO bi = myBxLo(myThid), myBxHi(myThid)
DO i=1,sNx
DO j=1,sNy
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_umask(i,j,bi,bj).eq.1.0) THEN
IF (streamice_hmask(i-1,j,bi,bj).eq.1.0.AND.
& streamice_hmask(i,j,bi,bj).eq.1.0) THEN
! geom_fac = sqrt(rA(i-1,j,bi,bj)*recip_rA(i,j,bi,bj)*
! & dxF(i,j,bi,bj)*recip_dxF(i-1,j,bi,bj))
taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) -
& 0.25 * dyG(i,j,bi,bj) *
& gravity *
& (H_streamice(i,j,bi,bj)*avg_density(i,j,bi,bj)+
& H_streamice(i-1,j,bi,bj)*avg_density(i-1,j,bi,bj)) *
& (surf_el_streamice(i,j,bi,bj)-
& surf_el_streamice(i-1,j,bi,bj))
CCC
taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) -
& streamice_density * gravity *
& streamice_bg_surf_slope_x * .25 * rA(i,j,bi,bj) *
& (H_streamice(i-1,j,bi,bj) + H_streamice(i,j,bi,bj))
CCC
ELSE IF (streamice_hmask(i-1,j,bi,bj).eq.1.0) THEN
IF (float_frac_streamice(i-1,j,bi,bj) .eq. 1.0) THEN
#ifdef USE_ALT_RLOW
IF (R_low_si(i-1,j,bi,bj) .lt. 0. _d 0) then
#else
IF (R_low(i-1,j,bi,bj) .lt. 0. _d 0) then
#endif
grd_below_sl = 1. _d 0
else
grd_below_sl = 0. _d 0
endif
taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) +
& 0.25 * dyG(i,j,bi,bj) *
& gravity *
& (avg_density(i-1,j,bi,bj) * H_streamice(i-1,j,bi,bj)**2-
#ifdef USE_ALT_RLOW
& streamice_density_ocean_avg * grd_below_sl *
& R_low_si(i-1,j,bi,bj)**2)
#else
& streamice_density_ocean_avg * grd_below_sl *
& R_low(i-1,j,bi,bj)**2)
#endif
ELSE
#ifdef STREAMICE_FIRN_CORRECTION
if (STREAMICE_apply_firn_correction) then
if (H_streamice(i-1,j,bi,bj).lt.firn_depth) then
taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) +
& 0.25 * dyG(i,j,bi,bj) *
& streamice_density_firn * gravity *
& (1-streamice_density_firn*i_rhow) *
& H_streamice(i-1,j,bi,bj)**2
else
taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) +
& 0.25 * dyG(i,j,bi,bj) * gravity * (
& streamice_density_firn * firn_depth**2 +
& (h_streamice(i-1,j,bi,bj)-firn_depth) *
& (streamice_density_firn*firn_depth+streamice_density*
& (h_streamice(i-1,j,bi,bj)-streamice_firn_correction)) -
& streamice_density**2*i_rhow*
& (h_streamice(i-1,j,bi,bj)-streamice_firn_correction)**2
& )
endif
else
#endif
taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) +
& 0.25 * dyG(i,j,bi,bj) *
& streamice_density * gravity *
& (1-streamice_density*i_rhow) *
& H_streamice(i-1,j,bi,bj)**2
#ifdef STREAMICE_FIRN_CORRECTION
endif
#endif
ENDIF
ELSE IF (streamice_hmask(i,j,bi,bj).eq.1.0) THEN
IF (float_frac_streamice(i,j,bi,bj) .eq. 1.0) THEN
#ifdef USE_ALT_RLOW
IF (R_low_si(i,j,bi,bj) .lt. 0. _d 0) then
#else
IF (R_low(i,j,bi,bj) .lt. 0. _d 0) then
#endif
grd_below_sl = 1. _d 0
else
grd_below_sl = 0. _d 0
endif
taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) -
& 0.25 * dyG(i,j,bi,bj) *
& gravity *
& (avg_density(i,j,bi,bj) * H_streamice(i,j,bi,bj)**2 -
#ifdef USE_ALT_RLOW
& streamice_density_ocean_avg * grd_below_sl *
& R_low_si(i,j,bi,bj)**2)
#else
& streamice_density_ocean_avg * grd_below_sl *
& R_low(i,j,bi,bj)**2)
#endif
ELSE
#ifdef STREAMICE_FIRN_CORRECTION
if (STREAMICE_apply_firn_correction) then
if (H_streamice(i,j,bi,bj).lt.firn_depth) then
taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) -
& 0.25 * dyG(i,j,bi,bj) *
& streamice_density_firn * gravity *
& (1-streamice_density_firn*i_rhow) *
& H_streamice(i,j,bi,bj)**2
else
taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) -
& 0.25 * dyG(i,j,bi,bj) * gravity * (
& streamice_density_firn * firn_depth**2 +
& (h_streamice(i,j,bi,bj)-firn_depth) *
& (streamice_density_firn*firn_depth+streamice_density*
& (h_streamice(i,j,bi,bj)-streamice_firn_correction)) -
& streamice_density**2*i_rhow*
& (h_streamice(i,j,bi,bj)-streamice_firn_correction)**2
& )
endif
else
#endif
taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) -
& 0.25 * dyG(i,j,bi,bj) *
& streamice_density * gravity *
& (1-streamice_density*i_rhow) *
& H_streamice(i,j,bi,bj)**2
#ifdef STREAMICE_FIRN_CORRECTION
endif
#endif
END
IF
END
IF
! cells below
IF (streamice_hmask(i-1,j-1,bi,bj).eq.1.0.AND.
& streamice_hmask(i,j-1,bi,bj).eq.1.0) THEN
! geom_fac = sqrt(rA(i-1,j-1,bi,bj)*recip_rA(i,j-1,bi,bj)*
! & dxF(i,j-1,bi,bj)*recip_dxF(i-1,j-1,bi,bj))
taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) -
& 0.25 * dyg(i,j-1,bi,bj) *
& gravity *
& (H_streamice(i,j-1,bi,bj)*avg_density(i,j-1,bi,bj)+
& H_streamice(i-1,j-1,bi,bj)*avg_density(i-1,j-1,bi,bj)) *
& (surf_el_streamice(i,j-1,bi,bj)-
& surf_el_streamice(i-1,j-1,bi,bj))
taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) -
& streamice_density * gravity *
& streamice_bg_surf_slope_x * .25 * rA(i,j-1,bi,bj) *
& (H_streamice(i-1,j-1,bi,bj) + H_streamice(i,j-1,bi,bj))
ELSE IF (streamice_hmask(i-1,j-1,bi,bj).eq.1.0) THEN
IF (float_frac_streamice(i-1,j-1,bi,bj) .eq. 1.0) THEN
#ifdef USE_ALT_RLOW
IF (R_low_si(i-1,j-1,bi,bj) .lt. 0. _d 0) then
#else
IF (R_low(i-1,j-1,bi,bj) .lt. 0. _d 0) then
#endif
grd_below_sl = 1. _d 0
else
grd_below_sl = 0. _d 0
endif
taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) +
& 0.25 * dyg(i,j-1,bi,bj) *
& gravity *
& (avg_density(i-1,j-1,bi,bj)*H_streamice(i-1,j-1,bi,bj)**2-
#ifdef USE_ALT_RLOW
& streamice_density_ocean_avg*grd_below_sl *
& R_low_si(i-1,j-1,bi,bj)**2)
#else
& streamice_density_ocean_avg * grd_below_sl *
& R_low(i-1,j-1,bi,bj)**2)
#endif
ELSE
#ifdef STREAMICE_FIRN_CORRECTION
if (STREAMICE_apply_firn_correction) then
if (H_streamice(i-1,j-1,bi,bj).lt.firn_depth) then
taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) +
& 0.25 * dyG(i,j-1,bi,bj) *
& streamice_density_firn * gravity *
& (1-streamice_density_firn*i_rhow) *
& H_streamice(i-1,j-1,bi,bj)**2
else
taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) +
& 0.25 * dyG(i,j-1,bi,bj) * gravity * (
& streamice_density_firn * firn_depth**2 +
& (h_streamice(i-1,j-1,bi,bj)-firn_depth) *
& (streamice_density_firn*firn_depth+streamice_density*
& (h_streamice(i-1,j-1,bi,bj)-streamice_firn_correction))-
& streamice_density**2*i_rhow*
& (h_streamice(i-1,j-1,bi,bj)-streamice_firn_correction)**2
& )
endif
else
#endif
taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) +
& 0.25 * dyG(i,j-1,bi,bj) *
& streamice_density * gravity *
& (1-streamice_density*i_rhow) *
& H_streamice(i-1,j-1,bi,bj)**2
#ifdef STREAMICE_FIRN_CORRECTION
endif
#endif
ENDIF
ELSE IF (streamice_hmask(i,j-1,bi,bj).eq.1.0) THEN
IF (float_frac_streamice(i,j-1,bi,bj) .eq. 1.0) THEN
#ifdef USE_ALT_RLOW
IF (R_low_si(i,j-1,bi,bj) .lt. 0. _d 0) then
#else
IF (R_low(i,j-1,bi,bj) .lt. 0. _d 0) then
#endif
grd_below_sl = 1. _d 0
else
grd_below_sl = 0. _d 0
endif
taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) -
& 0.25 * dyg(i,j-1,bi,bj) *
& gravity *
& (avg_density(i,j-1,bi,bj) * H_streamice(i,j-1,bi,bj)**2 -
#ifdef USE_ALT_RLOW
& streamice_density_ocean_avg * grd_below_sl *
& R_low_si(i,j-1,bi,bj)**2)
#else
& streamice_density_ocean_avg * grd_below_sl *
& R_low(i,j-1,bi,bj)**2)
#endif
ELSE
#ifdef STREAMICE_FIRN_CORRECTION
if (STREAMICE_apply_firn_correction) then
if (H_streamice(i,j-1,bi,bj).lt.firn_depth) then
taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) -
& 0.25 * dyG(i,j-1,bi,bj) *
& streamice_density_firn * gravity *
& (1-streamice_density_firn*i_rhow) *
& H_streamice(i,j-1,bi,bj)**2
else
taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) -
& 0.25 * dyG(i,j-1,bi,bj) * gravity * (
& streamice_density_firn * firn_depth**2 +
& (h_streamice(i,j-1,bi,bj)-firn_depth) *
& (streamice_density_firn*firn_depth+streamice_density*
& (h_streamice(i,j-1,bi,bj)-streamice_firn_correction))-
& streamice_density**2*i_rhow*
& (h_streamice(i,j-1,bi,bj)-streamice_firn_correction)**2
& )
endif
else
#endif
taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) -
& 0.25 * dyG(i,j-1,bi,bj) *
& streamice_density * gravity *
& (1-streamice_density*i_rhow) *
& H_streamice(i,j-1,bi,bj)**2
#ifdef STREAMICE_FIRN_CORRECTION
endif
#endif
END
IF
END
IF
END
IF ! if umask==1
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
IF (streamice_vmask(i,j,bi,bj).eq.1.0) THEN
IF (streamice_hmask(i,j-1,bi,bj).eq.1.0.AND.
& streamice_hmask(i,j,bi,bj).eq.1.0) THEN
! geom_fac = sqrt(rA(i,j-1,bi,bj)*recip_rA(i,j,bi,bj)*
! & dxF(i,j,bi,bj)*recip_dyF(i,j-1,bi,bj))
taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) -
& 0.25 * dxG(i,j,bi,bj) *
& gravity *
& (H_streamice(i,j,bi,bj)*avg_density(i,j,bi,bj)+
& H_streamice(i,j-1,bi,bj)*avg_density(i,j-1,bi,bj))*
& (surf_el_streamice(i,j,bi,bj)-
& surf_el_streamice(i,j-1,bi,bj))
taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) -
& streamice_density * gravity *
& streamice_bg_surf_slope_y * .25 * rA(i,j,bi,bj) *
& (H_streamice(i,j-1,bi,bj) + H_streamice(i,j,bi,bj))
ELSE IF (streamice_hmask(i,j-1,bi,bj).eq.1.0) THEN
IF (float_frac_streamice(i,j-1,bi,bj) .eq. 1.0) THEN
#ifdef USE_ALT_RLOW
IF (R_low_si(i,j-1,bi,bj) .lt. 0. _d 0) then
#else
IF (R_low(i,j-1,bi,bj) .lt. 0. _d 0) then
#endif
grd_below_sl = 1. _d 0
else
grd_below_sl = 0. _d 0
endif
taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) +
& 0.25 * dxG(i,j,bi,bj) *
& gravity *
& (avg_density(i,j-1,bi,bj) * H_streamice(i,j-1,bi,bj)**2 -
#ifdef USE_ALT_RLOW
& streamice_density_ocean_avg * grd_below_sl *
& R_low_si(i,j-1,bi,bj)**2)
#else
& streamice_density_ocean_avg * grd_below_sl *
& R_low(i,j-1,bi,bj)**2)
#endif
ELSE
#ifdef STREAMICE_FIRN_CORRECTION
if (STREAMICE_apply_firn_correction) then
if (H_streamice(i,j-1,bi,bj).lt.firn_depth) then
taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) +
& 0.25 * dxG(i,j,bi,bj) *
& streamice_density_firn * gravity *
& (1-streamice_density_firn*i_rhow) *
& H_streamice(i,j-1,bi,bj)**2
else
taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) +
& 0.25 * dxG(i,j,bi,bj) * gravity * (
& streamice_density_firn * firn_depth**2 +
& (h_streamice(i,j-1,bi,bj)-firn_depth) *
& (streamice_density_firn*firn_depth+streamice_density*
& (h_streamice(i,j-1,bi,bj)-streamice_firn_correction))-
& streamice_density**2*i_rhow*
& (h_streamice(i,j-1,bi,bj)-streamice_firn_correction)**2
& )
endif
else
#endif
taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) +
& 0.25 * dxG(i,j,bi,bj) *
& streamice_density * gravity *
& (1-streamice_density*i_rhow) *
& H_streamice(i,j-1,bi,bj)**2
#ifdef STREAMICE_FIRN_CORRECTION
endif
#endif
ENDIF
ELSE IF (streamice_hmask(i,j,bi,bj).eq.1.0) THEN
IF (float_frac_streamice(i,j,bi,bj) .eq. 1.0) THEN
#ifdef USE_ALT_RLOW
IF (R_low_si(i,j,bi,bj) .lt. 0. _d 0) then
#else
IF (R_low(i,j,bi,bj) .lt. 0. _d 0) then
#endif
grd_below_sl = 1. _d 0
else
grd_below_sl = 0. _d 0
endif
taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) -
& 0.25 * dxG(i,j,bi,bj) *
& gravity *
& (avg_density(i,j,bi,bj) * H_streamice(i,j,bi,bj)**2 -
#ifdef USE_ALT_RLOW
& streamice_density_ocean_avg * grd_below_sl *
& R_low_si(i,j,bi,bj)**2)
#else
& streamice_density_ocean_avg * grd_below_sl *
& R_low(i,j,bi,bj)**2)
#endif
ELSE
#ifdef STREAMICE_FIRN_CORRECTION
if (STREAMICE_apply_firn_correction) then
if (H_streamice(i,j,bi,bj).lt.firn_depth) then
taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) -
& 0.25 * dxG(i,j,bi,bj) *
& streamice_density_firn * gravity *
& (1-streamice_density_firn*i_rhow) *
& H_streamice(i,j,bi,bj)**2
else
taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) -
& 0.25 * dxG(i,j,bi,bj) * gravity * (
& streamice_density_firn * firn_depth**2 +
& (h_streamice(i,j,bi,bj)-firn_depth) *
& (streamice_density_firn*firn_depth+streamice_density*
& (h_streamice(i,j,bi,bj)-streamice_firn_correction))-
& streamice_density**2*i_rhow*
& (h_streamice(i,j,bi,bj)-streamice_firn_correction)**2
& )
endif
else
#endif
taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) -
& 0.25 * dxG(i,j,bi,bj) *
& streamice_density * gravity *
& (1-streamice_density*i_rhow) *
& H_streamice(i,j,bi,bj)**2
#ifdef STREAMICE_FIRN_CORRECTION
endif
#endif
END
IF
END
IF
! cells to left
IF (streamice_hmask(i-1,j-1,bi,bj).eq.1.0.AND.
& streamice_hmask(i-1,j,bi,bj).eq.1.0) THEN
! geom_fac = sqrt(rA(i-1,j-1,bi,bj)*recip_rA(i-1,j,bi,bj)*
! & dxF(i-1,j,bi,bj)*recip_dxF(i-1,j-1,bi,bj))
taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) -
& 0.25 * dxG(i-1,j,bi,bj) *
& gravity *
& (H_streamice(i-1,j,bi,bj)*avg_density(i-1,j,bi,bj)+
& H_streamice(i-1,j-1,bi,bj)*avg_density(i-1,j-1,bi,bj))*
& (surf_el_streamice(i-1,j,bi,bj)-
& surf_el_streamice(i-1,j-1,bi,bj))
taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) -
& streamice_density * gravity *
& streamice_bg_surf_slope_y * .25 * rA(i-1,j,bi,bj) *
& (H_streamice(i-1,j-1,bi,bj) + H_streamice(i-1,j,bi,bj))
ELSE IF (streamice_hmask(i-1,j-1,bi,bj).eq.1.0) THEN
IF (float_frac_streamice(i-1,j-1,bi,bj) .eq. 1.0) THEN
#ifdef USE_ALT_RLOW
IF (R_low_si(i-1,j-1,bi,bj) .lt. 0. _d 0) then
#else
IF (R_low(i-1,j-1,bi,bj) .lt. 0. _d 0) then
#endif
grd_below_sl = 1. _d 0
else
grd_below_sl = 0. _d 0
endif
taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) +
& 0.25 * dxG(i-1,j,bi,bj) *
& gravity *
& (avg_density(i-1,j-1,bi,bj)*H_streamice(i-1,j-1,bi,bj)**2-
#ifdef USE_ALT_RLOW
& streamice_density_ocean_avg*grd_below_sl *
& R_low_si(i-1,j-1,bi,bj)**2)
#else
& streamice_density_ocean_avg * grd_below_sl *
& R_low(i-1,j-1,bi,bj)**2)
#endif
ELSE
#ifdef STREAMICE_FIRN_CORRECTION
if (STREAMICE_apply_firn_correction) then
if (H_streamice(i-1,j-1,bi,bj).lt.firn_depth) then
taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) +
& 0.25 * dxG(i-1,j,bi,bj) *
& streamice_density_firn * gravity *
& (1-streamice_density_firn*i_rhow) *
& H_streamice(i-1,j-1,bi,bj)**2
else
taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) +
& 0.25 * dxG(i-1,j,bi,bj) * gravity * (
& streamice_density_firn * firn_depth**2 +
& (h_streamice(i-1,j-1,bi,bj)-firn_depth) *
& (streamice_density_firn*firn_depth+streamice_density*
& (h_streamice(i-1,j-1,bi,bj)-streamice_firn_correction))-
& streamice_density**2*i_rhow*
& (h_streamice(i-1,j-1,bi,bj)-streamice_firn_correction)**2
& )
endif
else
#endif
taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) +
& 0.25 * dxG(i-1,j,bi,bj) *
& streamice_density * gravity *
& (1-streamice_density*i_rhow) *
& H_streamice(i-1,j-1,bi,bj)**2
#ifdef STREAMICE_FIRN_CORRECTION
endif
#endif
ENDIF
ELSE IF (streamice_hmask(i-1,j,bi,bj).eq.1.0) THEN
IF (float_frac_streamice(i-1,j,bi,bj) .eq. 1.0) THEN
#ifdef USE_ALT_RLOW
IF (R_low_si(i-1,j,bi,bj) .lt. 0. _d 0) then
#else
IF (R_low(i-1,j,bi,bj) .lt. 0. _d 0) then
#endif
grd_below_sl = 1. _d 0
else
grd_below_sl = 0. _d 0
endif
taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) -
& 0.25 * dxG(i-1,j,bi,bj) *
& gravity *
& (avg_density(i-1,j,bi,bj) * H_streamice(i-1,j,bi,bj)**2 -
#ifdef USE_ALT_RLOW
& streamice_density_ocean_avg * grd_below_sl *
& R_low_si(i-1,j,bi,bj)**2)
#else
& streamice_density_ocean_avg * grd_below_sl *
& R_low(i-1,j,bi,bj)**2)
#endif
ELSE
#ifdef STREAMICE_FIRN_CORRECTION
if (STREAMICE_apply_firn_correction) then
if (H_streamice(i-1,j,bi,bj).lt.firn_depth) then
taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) -
& 0.25 * dxG(i-1,j,bi,bj) *
& streamice_density_firn * gravity *
& (1-streamice_density_firn*i_rhow) *
& H_streamice(i-1,j,bi,bj)**2
else
taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) -
& 0.25 * dxG(i-1,j,bi,bj) * gravity * (
& streamice_density_firn * firn_depth**2 +
& (h_streamice(i-1,j,bi,bj)-firn_depth) *
& (streamice_density_firn*firn_depth+streamice_density*
& (h_streamice(i-1,j,bi,bj)-streamice_firn_correction))-
& streamice_density**2*i_rhow*
& (h_streamice(i-1,j,bi,bj)-streamice_firn_correction)**2
& )
endif
else
#endif
taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) -
& 0.25 * dxG(i-1,j,bi,bj) *
& streamice_density * gravity *
& (1-streamice_density*i_rhow) *
& H_streamice(i-1,j,bi,bj)**2
#ifdef STREAMICE_FIRN_CORRECTION
endif
#endif
END
IF
END
IF
END
IF ! if vmask ==1
ENDDO
ENDDO
ENDDO
ENDDO
! taudx_SI (1,1,1,1) = taudx_SI (1,1,1,1) +
! & streamice_v_normal_pert (1,1,1,1)
! call write_fld_xy_rl("TAUDX_SI","",taudx_si,0,mythid)
! call write_fld_xy_rl("TAUDY_SI","",taudy_si,0,mythid)
#endif
RETURN
END