Skip to content

Commit

Permalink
ice_dyn_vp: remove unused subroutines
Browse files Browse the repository at this point in the history 
Remove subroutines 'stress_prime_vpOLD', 'matvecOLD' and 'precond_diag'.

'stress_prime_vpOLD' and 'matvecOLD' have been unused since the creation
of the 'matvec' subroutine.

'precond_diag' has been unused since we deactivated the legacy FGMRES
solver.
  • Loading branch information
@phil-blain
phil-blain committed Jul 15, 2020
1 parent adb62f0 commit 083b611
Showing 1 changed file with 0 additions and 355 deletions.
355 changes: 0 additions & 355 deletions cicecore/cicedynB/dynamics/ice_dyn_vp.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1292,237 +1292,6 @@ end subroutine calc_zeta_Pr

!=======================================================================

! Computes VP stress without the rep. pressure Pr (included in b vector)

subroutine stress_prime_vpOLD (nx_block, ny_block, &
icellt, &
indxti, indxtj, &
uvel, vvel, &
dxt, dyt, &
dxhy, dyhx, &
cxp, cyp, &
cxm, cym, &
zetaD, &
str )

use ice_dyn_shared, only: strain_rates

integer (kind=int_kind), intent(in) :: &
nx_block, ny_block, & ! block dimensions
icellt ! no. of cells where icetmask = 1

integer (kind=int_kind), dimension (nx_block*ny_block), &
intent(in) :: &
indxti , & ! compressed index in i-direction
indxtj ! compressed index in j-direction

real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
uvel , & ! x-component of velocity (m/s)
vvel , & ! y-component of velocity (m/s)
dxt , & ! width of T-cell through the middle (m)
dyt , & ! height of T-cell through the middle (m)
dxhy , & ! 0.5*(HTE - HTE)
dyhx , & ! 0.5*(HTN - HTN)
cyp , & ! 1.5*HTE - 0.5*HTE
cxp , & ! 1.5*HTN - 0.5*HTN
cym , & ! 0.5*HTE - 1.5*HTE
cxm ! 0.5*HTN - 1.5*HTN

real (kind=dbl_kind), dimension(nx_block,ny_block,4), &
intent(in) :: &
zetaD ! 2*zeta

real (kind=dbl_kind), dimension(nx_block,ny_block,8), &
intent(out) :: &
str ! stress combinations

! local variables

integer (kind=int_kind) :: &
i, j, ij

real (kind=dbl_kind) :: &
divune, divunw, divuse, divusw , & ! divergence
tensionne, tensionnw, tensionse, tensionsw, & ! tension
shearne, shearnw, shearse, shearsw , & ! shearing
Deltane, Deltanw, Deltase, Deltasw , & ! Delt
ssigpn, ssigps, ssigpe, ssigpw , &
ssigmn, ssigms, ssigme, ssigmw , &
ssig12n, ssig12s, ssig12e, ssig12w , &
ssigp1, ssigp2, ssigm1, ssigm2, ssig121, ssig122, &
csigpne, csigpnw, csigpse, csigpsw , &
csigmne, csigmnw, csigmse, csigmsw , &
csig12ne, csig12nw, csig12se, csig12sw , &
str12ew, str12we, str12ns, str12sn , &
strp_tmp, strm_tmp

real (kind=dbl_kind) :: &
stressp_1, stressp_2, stressp_3, stressp_4 , & ! sigma11+sigma22 (without Pr)
stressm_1, stressm_2, stressm_3, stressm_4 , & ! sigma11-sigma22
stress12_1,stress12_2,stress12_3,stress12_4 ! sigma12

character(len=*), parameter :: subname = '(stress_prime_vpOLD)'

!-----------------------------------------------------------------
! Initialize
!-----------------------------------------------------------------

str(:,:,:) = c0

!DIR$ CONCURRENT !Cray
!cdir nodep !NEC
!ocl novrec !Fujitsu

do ij = 1, icellt
i = indxti(ij)
j = indxtj(ij)

!-----------------------------------------------------------------
! strain rates
! NOTE these are actually strain rates * area (m^2/s)
!-----------------------------------------------------------------
call strain_rates (nx_block, ny_block, &
i, j, &
uvel, vvel, &
dxt, dyt, &
cxp, cyp, &
cxm, cym, &
divune, divunw, &
divuse, divusw, &
tensionne, tensionnw, &
tensionse, tensionsw, &
shearne, shearnw, &
shearse, shearsw, &
Deltane, Deltanw, &
Deltase, Deltasw )

!-----------------------------------------------------------------
! the stresses ! kg/s^2
! (1) northeast, (2) northwest, (3) southwest, (4) southeast
! JFL commented part of stressp is for the rep pressure Pr
!-----------------------------------------------------------------

stressp_1 = zetaD(i,j,1)*(divune*(c1+Ktens))! - Deltane*(c1-Ktens))
stressp_2 = zetaD(i,j,2)*(divunw*(c1+Ktens))! - Deltanw*(c1-Ktens))
stressp_3 = zetaD(i,j,3)*(divusw*(c1+Ktens))! - Deltasw*(c1-Ktens))
stressp_4 = zetaD(i,j,4)*(divuse*(c1+Ktens))! - Deltase*(c1-Ktens))

stressm_1 = zetaD(i,j,1)*tensionne*(c1+Ktens)*ecci
stressm_2 = zetaD(i,j,2)*tensionnw*(c1+Ktens)*ecci
stressm_3 = zetaD(i,j,3)*tensionsw*(c1+Ktens)*ecci
stressm_4 = zetaD(i,j,4)*tensionse*(c1+Ktens)*ecci

stress12_1 = zetaD(i,j,1)*shearne*p5*(c1+Ktens)*ecci
stress12_2 = zetaD(i,j,2)*shearnw*p5*(c1+Ktens)*ecci
stress12_3 = zetaD(i,j,3)*shearsw*p5*(c1+Ktens)*ecci
stress12_4 = zetaD(i,j,4)*shearse*p5*(c1+Ktens)*ecci

!-----------------------------------------------------------------
! combinations of the stresses for the momentum equation ! kg/s^2
!-----------------------------------------------------------------

ssigpn = stressp_1 + stressp_2
ssigps = stressp_3 + stressp_4
ssigpe = stressp_1 + stressp_4
ssigpw = stressp_2 + stressp_3
ssigp1 =(stressp_1 + stressp_3)*p055
ssigp2 =(stressp_2 + stressp_4)*p055

ssigmn = stressm_1 + stressm_2
ssigms = stressm_3 + stressm_4
ssigme = stressm_1 + stressm_4
ssigmw = stressm_2 + stressm_3
ssigm1 =(stressm_1 + stressm_3)*p055
ssigm2 =(stressm_2 + stressm_4)*p055

ssig12n = stress12_1 + stress12_2
ssig12s = stress12_3 + stress12_4
ssig12e = stress12_1 + stress12_4
ssig12w = stress12_2 + stress12_3
ssig121 =(stress12_1 + stress12_3)*p111
ssig122 =(stress12_2 + stress12_4)*p111

csigpne = p111*stressp_1 + ssigp2 + p027*stressp_3
csigpnw = p111*stressp_2 + ssigp1 + p027*stressp_4
csigpsw = p111*stressp_3 + ssigp2 + p027*stressp_1
csigpse = p111*stressp_4 + ssigp1 + p027*stressp_2

csigmne = p111*stressm_1 + ssigm2 + p027*stressm_3
csigmnw = p111*stressm_2 + ssigm1 + p027*stressm_4
csigmsw = p111*stressm_3 + ssigm2 + p027*stressm_1
csigmse = p111*stressm_4 + ssigm1 + p027*stressm_2

csig12ne = p222*stress12_1 + ssig122 &
+ p055*stress12_3
csig12nw = p222*stress12_2 + ssig121 &
+ p055*stress12_4
csig12sw = p222*stress12_3 + ssig122 &
+ p055*stress12_1
csig12se = p222*stress12_4 + ssig121 &
+ p055*stress12_2

str12ew = p5*dxt(i,j)*(p333*ssig12e + p166*ssig12w)
str12we = p5*dxt(i,j)*(p333*ssig12w + p166*ssig12e)
str12ns = p5*dyt(i,j)*(p333*ssig12n + p166*ssig12s)
str12sn = p5*dyt(i,j)*(p333*ssig12s + p166*ssig12n)

!-----------------------------------------------------------------
! for dF/dx (u momentum)
!-----------------------------------------------------------------
strp_tmp = p25*dyt(i,j)*(p333*ssigpn + p166*ssigps)
strm_tmp = p25*dyt(i,j)*(p333*ssigmn + p166*ssigms)

! northeast (i,j)
str(i,j,1) = -strp_tmp - strm_tmp - str12ew &
+ dxhy(i,j)*(-csigpne + csigmne) + dyhx(i,j)*csig12ne

! northwest (i+1,j)
str(i,j,2) = strp_tmp + strm_tmp - str12we &
+ dxhy(i,j)*(-csigpnw + csigmnw) + dyhx(i,j)*csig12nw

strp_tmp = p25*dyt(i,j)*(p333*ssigps + p166*ssigpn)
strm_tmp = p25*dyt(i,j)*(p333*ssigms + p166*ssigmn)

! southeast (i,j+1)
str(i,j,3) = -strp_tmp - strm_tmp + str12ew &
+ dxhy(i,j)*(-csigpse + csigmse) + dyhx(i,j)*csig12se

! southwest (i+1,j+1)
str(i,j,4) = strp_tmp + strm_tmp + str12we &
+ dxhy(i,j)*(-csigpsw + csigmsw) + dyhx(i,j)*csig12sw

!-----------------------------------------------------------------
! for dF/dy (v momentum)
!-----------------------------------------------------------------
strp_tmp = p25*dxt(i,j)*(p333*ssigpe + p166*ssigpw)
strm_tmp = p25*dxt(i,j)*(p333*ssigme + p166*ssigmw)

! northeast (i,j)
str(i,j,5) = -strp_tmp + strm_tmp - str12ns &
- dyhx(i,j)*(csigpne + csigmne) + dxhy(i,j)*csig12ne

! southeast (i,j+1)
str(i,j,6) = strp_tmp - strm_tmp - str12sn &
- dyhx(i,j)*(csigpse + csigmse) + dxhy(i,j)*csig12se

strp_tmp = p25*dxt(i,j)*(p333*ssigpw + p166*ssigpe)
strm_tmp = p25*dxt(i,j)*(p333*ssigmw + p166*ssigme)

! northwest (i+1,j)
str(i,j,7) = -strp_tmp + strm_tmp + str12ns &
- dyhx(i,j)*(csigpnw + csigmnw) + dxhy(i,j)*csig12nw

! southwest (i+1,j+1)
str(i,j,8) = strp_tmp - strm_tmp + str12sn &
- dyhx(i,j)*(csigpsw + csigmsw) + dxhy(i,j)*csig12sw

enddo ! ij

end subroutine stress_prime_vpOLD


!=======================================================================

! Computes the VP stress (as diagnostic)

subroutine stress_vp (nx_block, ny_block, &
Expand Down Expand Up @@ -1753,92 +1522,6 @@ end subroutine calc_seabed_stress

!=======================================================================

! OLD matrix vector product A(u,v) * (u,v)
! Au = A(u,v)_[x] * uvel (x components of A(u,v) * (u,v))
! Av = A(u,v)_[y] * vvel (y components of A(u,v) * (u,v))

subroutine matvecOLD (nx_block, ny_block, &
icellu, &
indxui, indxuj, &
str, &
vrel, &
umassdti, fm, &
uarear, Cb, &
uvel, vvel, &
Au, Av)

integer (kind=int_kind), intent(in) :: &
nx_block, ny_block, & ! block dimensions
icellu ! total count when iceumask is true

integer (kind=int_kind), dimension (nx_block*ny_block), &
intent(in) :: &
indxui , & ! compressed index in i-direction
indxuj ! compressed index in j-direction

real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
vrel, & ! coefficient for tauw
Cb, & ! coefficient for basal stress
umassdti, & ! mass of U-cell/dt (kg/m^2 s)
fm , & ! Coriolis param. * mass in U-cell (kg/s)
uarear ! 1/uarea

real (kind=dbl_kind), dimension(nx_block,ny_block,8), &
intent(in) :: &
str

real (kind=dbl_kind), dimension (nx_block,ny_block), &
intent(in) :: &
uvel , & ! x-component of velocity (m/s)
vvel ! y-component of velocity (m/s)

real (kind=dbl_kind), dimension (nx_block,ny_block), &
intent(inout) :: &
Au , & ! matvec, Fx = bx - Au (N/m^2)
Av ! matvec, Fy = by - Av (N/m^2)

! local variables

integer (kind=int_kind) :: &
i, j, ij

real (kind=dbl_kind) :: &
utp, vtp , & ! utp = uvel, vtp = vvel
ccaimp,ccb , & ! intermediate variables
strintx, strinty

character(len=*), parameter :: subname = '(matvecOLD)'

!-----------------------------------------------------------------
! integrate the momentum equation
!-----------------------------------------------------------------

do ij =1, icellu
i = indxui(ij)
j = indxuj(ij)

utp = uvel(i,j)
vtp = vvel(i,j)

ccaimp = umassdti(i,j) + vrel(i,j) * cosw + Cb(i,j) ! kg/m^2 s

ccb = fm(i,j) + sign(c1,fm(i,j)) * vrel(i,j) * sinw ! kg/m^2 s

! divergence of the internal stress tensor
strintx = uarear(i,j)* &
(str(i,j,1) + str(i+1,j,2) + str(i,j+1,3) + str(i+1,j+1,4))
strinty = uarear(i,j)* &
(str(i,j,5) + str(i,j+1,6) + str(i+1,j,7) + str(i+1,j+1,8))

Au(i,j) = ccaimp*utp - ccb*vtp - strintx
Av(i,j) = ccaimp*vtp + ccb*utp - strinty

enddo ! ij

end subroutine matvecOLD

!=======================================================================

! Computes the matrix vector product A(u,v) * (u,v)
! Au = A(u,v)_[x] * uvel (x components of A(u,v) * (u,v))
! Av = A(u,v)_[y] * vvel (y components of A(u,v) * (u,v))
Expand Down Expand Up @@ -2803,44 +2486,6 @@ end subroutine formDiag_step2

!=======================================================================

! Diagonal (Jacobi) preconditioner for the legacy FGMRES driver

subroutine precond_diag (ntot, &
diagvec, &
wk1, wk2)

integer (kind=int_kind), intent(in) :: &
ntot ! size of problem for fgmres

real (kind=dbl_kind), dimension (ntot), intent(in) :: &
diagvec, wk1

real (kind=dbl_kind), dimension (ntot), intent(out) :: &
wk2

! local variables

integer (kind=int_kind) :: &
i

character(len=*), parameter :: subname = '(precond_diag)'

!-----------------------------------------------------------------
! form vector (converts from max_blocks arrays to single vector
!-----------------------------------------------------------------

wk2(:)=c0

do i=1, ntot

wk2(i) = wk1(i)/diagvec(i)

enddo! i

end subroutine precond_diag

!=======================================================================

! Compute squared l^2 norm of a grid function (tpu,tpv)

subroutine calc_L2norm_squared (nx_block, ny_block, &
Expand Down

0 comments on commit 083b611

Please sign in to comment.