diff --git a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 index a2ce00f59..76163b867 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 @@ -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, & @@ -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)) @@ -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, &