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