From 4147db49dbc93ba93cbb5560b1c299f96a98f4d7 Mon Sep 17 00:00:00 2001 From: JFLemieux73 <31927797+JFLemieux73@users.noreply.github.com> Date: Thu, 21 Oct 2021 11:14:25 -0400 Subject: [PATCH] ice_dyn_vp: express rheology in terms of bulk and shear viscosities (#647) * Added array for eta and changed name of zetaD...BFB * Added array for replacement pressure...BFB * Replacement pressure now used on calc...BFB * eta is now used in the calc...BFB * Ktens included in zeta,eta and rep_prs...roundoff errors...BFB only if Ktens=0 * Small modif to viscous coeff subroutine for passing logical capping * Modifs to viscous_coeffs suroutine to be used for VP and EVP * Further modifs to viscous_coeff subroutine * cosmetic change: order of calc is ne-nw-sw-se * vp solver also uses viscous_coeffs_and_rep_pressure subroutine * Minor modifications following PR review --- cicecore/cicedynB/dynamics/ice_dyn_evp.F90 | 12 +- cicecore/cicedynB/dynamics/ice_dyn_shared.F90 | 57 +++-- cicecore/cicedynB/dynamics/ice_dyn_vp.F90 | 242 ++++++++++-------- 3 files changed, 175 insertions(+), 136 deletions(-) diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index 775f2caf1..8f3fc4910 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 @@ -656,6 +656,8 @@ subroutine stress (nx_block, ny_block, & str12ew, str12we, str12ns, str12sn , & strp_tmp, strm_tmp, tmp + logical :: capping ! of the viscous coef + character(len=*), parameter :: subname = '(stress)' !----------------------------------------------------------------- @@ -663,6 +665,7 @@ subroutine stress (nx_block, ny_block, & !----------------------------------------------------------------- str(:,:,:) = c0 + capping = .true. ! could be later included in ice_in do ij = 1, icellt i = indxti(ij) @@ -694,13 +697,14 @@ subroutine stress (nx_block, ny_block, & call viscous_coeffs_and_rep_pressure (strength(i,j), tinyarea(i,j),& Deltane, Deltanw, & - Deltase, Deltasw, & + Deltasw, Deltase, & zetax2ne, zetax2nw, & - zetax2se, zetax2sw, & + zetax2sw, zetax2se, & etax2ne, etax2nw, & - etax2se, etax2sw, & + etax2sw, etax2se, & rep_prsne, rep_prsnw, & - rep_prsse, rep_prssw ) + rep_prssw, rep_prsse, & + capping) !----------------------------------------------------------------- ! the stresses ! kg/s^2 diff --git a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 index 167bc6a2c..5e14d9686 100755 --- a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 @@ -1375,53 +1375,64 @@ end subroutine strain_rates subroutine viscous_coeffs_and_rep_pressure (strength, tinyarea, & Deltane, Deltanw, & - Deltase, Deltasw, & + Deltasw, Deltase, & zetax2ne, zetax2nw, & - zetax2se, zetax2sw, & + zetax2sw, zetax2se, & etax2ne, etax2nw, & - etax2se, etax2sw, & + etax2sw, etax2se, & rep_prsne, rep_prsnw,& - rep_prsse, rep_prssw ) + rep_prssw, rep_prsse,& + capping) real (kind=dbl_kind), intent(in):: & strength, tinyarea ! at the t-point real (kind=dbl_kind), intent(in):: & - Deltane, Deltanw, Deltase, Deltasw ! Delta at each corner + Deltane, Deltanw, Deltasw, Deltase ! Delta at each corner + logical , intent(in):: capping + real (kind=dbl_kind), intent(out):: & - zetax2ne, zetax2nw, zetax2se, zetax2sw, & ! zetax2 at each corner - etax2ne, etax2nw, etax2se, etax2sw, & ! etax2 at each corner - rep_prsne, rep_prsnw, rep_prsse, rep_prssw ! replacement pressure + zetax2ne, zetax2nw, zetax2sw, zetax2se, & ! zetax2 at each corner + etax2ne, etax2nw, etax2sw, etax2se, & ! etax2 at each corner + rep_prsne, rep_prsnw, rep_prssw, rep_prsse ! replacement pressure ! local variables real (kind=dbl_kind) :: & - tmpcalc + tmpcalcne, tmpcalcnw, tmpcalcsw, tmpcalcse ! NOTE: for comp. efficiency 2 x zeta and 2 x eta are used in the code ! if (trim(yield_curve) == 'ellipse') then - tmpcalc = strength/max(Deltane,tinyarea) ! northeast - zetax2ne = (c1+Ktens)*tmpcalc - rep_prsne = (c1-Ktens)*tmpcalc*Deltane + if (capping) then + tmpcalcne = strength/max(Deltane,tinyarea) + tmpcalcnw = strength/max(Deltanw,tinyarea) + tmpcalcsw = strength/max(Deltasw,tinyarea) + tmpcalcse = strength/max(Deltase,tinyarea) + else + tmpcalcne = strength/(Deltane + tinyarea) + tmpcalcnw = strength/(Deltanw + tinyarea) + tmpcalcsw = strength/(Deltasw + tinyarea) + tmpcalcse = strength/(Deltase + tinyarea) + endif + + zetax2ne = (c1+Ktens)*tmpcalcne ! northeast + rep_prsne = (c1-Ktens)*tmpcalcne*Deltane etax2ne = ecci*zetax2ne ! CHANGE FOR e_plasticpot - tmpcalc = strength/max(Deltanw,tinyarea) ! northwest - zetax2nw = (c1+Ktens)*tmpcalc - rep_prsnw = (c1-Ktens)*tmpcalc*Deltanw + zetax2nw = (c1+Ktens)*tmpcalcnw ! northwest + rep_prsnw = (c1-Ktens)*tmpcalcnw*Deltanw etax2nw = ecci*zetax2nw ! CHANGE FOR e_plasticpot - tmpcalc = strength/max(Deltase,tinyarea) ! southeast - zetax2se = (c1+Ktens)*tmpcalc - rep_prsse = (c1-Ktens)*tmpcalc*Deltase - etax2se = ecci*zetax2se ! CHANGE FOR e_plasticpot - - tmpcalc = strength/max(Deltasw,tinyarea) ! southwest - zetax2sw = (c1+Ktens)*tmpcalc - rep_prssw = (c1-Ktens)*tmpcalc*Deltasw + zetax2sw = (c1+Ktens)*tmpcalcsw ! southwest + rep_prssw = (c1-Ktens)*tmpcalcsw*Deltasw etax2sw = ecci*zetax2sw ! CHANGE FOR e_plasticpot + zetax2se = (c1+Ktens)*tmpcalcse ! southeast + rep_prsse = (c1-Ktens)*tmpcalcse*Deltase + etax2se = ecci*zetax2se ! CHANGE FOR e_plasticpot + ! else ! endif diff --git a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 index 860865dba..1a6c68548 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 @@ -234,8 +234,10 @@ subroutine implicit_solver (dt) umassdti ! mass of U-cell/dte (kg/m^2 s) real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks,4):: & - zetaD ! zetaD = 2zeta (viscous coeff) - + zetax2 , & ! zetax2 = 2*zeta (bulk viscous coeff) + etax2 , & ! etax2 = 2*eta (shear viscous coeff) + rep_prs ! replacement pressure + logical (kind=log_kind) :: calc_strair integer (kind=int_kind), dimension (nx_block,ny_block,max_blocks) :: & @@ -488,8 +490,9 @@ subroutine implicit_solver (dt) bxfix , byfix , & umassdti, sol , & fpresx , fpresy, & - zetaD , Cb , & - halo_info_mask) + zetax2 , etax2 , & + rep_prs , & + Cb, halo_info_mask) !----------------------------------------------------------------- ! End of nonlinear iteration !----------------------------------------------------------------- @@ -510,7 +513,8 @@ subroutine implicit_solver (dt) dxt (:,:,iblk), dyt (:,:,iblk), & cxp (:,:,iblk), cyp (:,:,iblk), & cxm (:,:,iblk), cym (:,:,iblk), & - zetaD (:,:,iblk,:), & + zetax2 (:,:,iblk,:), etax2 (:,:,iblk,:), & + rep_prs (:,:,iblk,:), & stressp_1 (:,:,iblk), stressp_2 (:,:,iblk), & stressp_3 (:,:,iblk), stressp_4 (:,:,iblk), & stressm_1 (:,:,iblk), stressm_2 (:,:,iblk), & @@ -671,8 +675,9 @@ subroutine anderson_solver (icellt , icellu, & bxfix , byfix , & umassdti, sol , & fpresx , fpresy, & - zetaD , Cb , & - halo_info_mask) + zetax2 , etax2 , & + rep_prs , & + Cb, halo_info_mask) use ice_arrays_column, only: Cdn_ocn use ice_blocks, only: nx_block, ny_block @@ -708,8 +713,10 @@ subroutine anderson_solver (icellt , icellu, & umassdti ! mass of U-cell/dte (kg/m^2 s) real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks,4), intent(out) :: & - zetaD ! zetaD = 2zeta (viscous coeff) - + zetax2 , & ! zetax2 = 2*zeta (bulk viscous coeff) + etax2 , & ! etax2 = 2*eta (shear viscous coeff) + rep_prs ! replacement pressure + type (ice_halo), intent(in) :: & halo_info_mask ! ghost cell update info for masked halo @@ -805,7 +812,7 @@ subroutine anderson_solver (icellt , icellu, & do it_nl = 0, maxits_nonlin ! nonlinear iteration loop ! Compute quantities needed for computing PDE residual and g(x) (fixed point map) !----------------------------------------------------------------- - ! Calc zetaD, dPr/dx, dPr/dy, Cb and vrel = f(uprev_k, vprev_k) + ! Calc zetax2, etax2, dPr/dx, dPr/dy, Cb and vrel = f(uprev_k, vprev_k) !----------------------------------------------------------------- !$OMP PARALLEL DO PRIVATE(iblk) do iblk = 1, nblocks @@ -828,9 +835,9 @@ subroutine anderson_solver (icellt , icellu, & dxhy (:,:,iblk), dyhx (:,:,iblk), & cxp (:,:,iblk), cyp (:,:,iblk), & cxm (:,:,iblk), cym (:,:,iblk), & - tinyarea (:,:,iblk), & - strength (:,:,iblk), zetaD (:,:,iblk,:), & - stress_Pr (:,:,:)) + tinyarea (:,:,iblk), strength (:,:,iblk),& + zetax2 (:,:,iblk,:), etax2 (:,:,iblk,:),& + rep_prs(:,:,iblk,:), stress_Pr (:,:,:)) call calc_vrel_Cb (nx_block , ny_block , & icellu (iblk), Cdn_ocn (:,:,iblk), & @@ -861,7 +868,7 @@ subroutine anderson_solver (icellt , icellu, & cxm (:,:,iblk) , cym (:,:,iblk), & uprev_k (:,:,iblk) , vprev_k (:,:,iblk), & vrel (:,:,iblk) , Cb (:,:,iblk), & - zetaD (:,:,iblk,:), & + zetax2 (:,:,iblk,:), etax2 (:,:,iblk,:), & umassdti (:,:,iblk) , fm (:,:,iblk), & uarear (:,:,iblk) , & Au (:,:,iblk) , Av (:,:,iblk)) @@ -908,14 +915,15 @@ subroutine anderson_solver (icellt , icellu, & !$OMP PARALLEL DO PRIVATE(iblk) do iblk = 1, nblocks ! first compute diagonal contributions due to rheology term - call formDiag_step1 (nx_block , ny_block , & - icellu (iblk) , & - indxui (:,iblk) , indxuj(:,iblk), & - dxt (:,:,iblk) , dyt (:,:,iblk), & - dxhy (:,:,iblk) , dyhx(:,:,iblk), & - cxp (:,:,iblk) , cyp (:,:,iblk), & - cxm (:,:,iblk) , cym (:,:,iblk), & - zetaD (:,:,iblk,:), diag_rheo(:,:,:)) + call formDiag_step1 (nx_block , ny_block , & + icellu (iblk) , & + indxui (:,iblk) , indxuj(:,iblk), & + dxt (:,:,iblk) , dyt (:,:,iblk), & + dxhy (:,:,iblk) , dyhx(:,:,iblk), & + cxp (:,:,iblk) , cyp (:,:,iblk), & + cxm (:,:,iblk) , cym (:,:,iblk), & + zetax2 (:,:,iblk,:), etax2(:,:,iblk,:), & + diag_rheo(:,:,:)) ! second compute the full diagonal call formDiag_step2 (nx_block , ny_block , & icellu (iblk), & @@ -929,7 +937,7 @@ subroutine anderson_solver (icellt , icellu, & endif ! FGMRES linear solver - call fgmres (zetaD , & + call fgmres (zetax2 , etax2 , & Cb , vrel , & umassdti , & halo_info_mask, & @@ -1119,7 +1127,7 @@ end subroutine anderson_solver !======================================================================= -! Computes the viscous coefficients (in fact zetaD=2*zeta) and dPr/dx, dPr/dy +! Computes the viscous coefficients and dPr/dx, dPr/dy subroutine calc_zeta_dPr (nx_block, ny_block, & icellt , & @@ -1129,11 +1137,11 @@ subroutine calc_zeta_dPr (nx_block, ny_block, & dxhy , dyhx , & cxp , cyp , & cxm , cym , & - tinyarea, & - strength, zetaD , & - stPr) + tinyarea, strength, & + zetax2 , etax2 , & + rep_prs , stPr) - use ice_dyn_shared, only: strain_rates + use ice_dyn_shared, only: strain_rates, viscous_coeffs_and_rep_pressure integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions @@ -1158,8 +1166,10 @@ subroutine calc_zeta_dPr (nx_block, ny_block, & tinyarea ! min_strain_rate*tarea real (kind=dbl_kind), dimension(nx_block,ny_block,4), intent(out) :: & - zetaD ! 2*zeta - + zetax2 , & ! zetax2 = 2*zeta (bulk viscous coeff) + etax2 , & ! etax2 = 2*eta (shear viscous coeff) + rep_prs ! replacement pressure + real (kind=dbl_kind), dimension(nx_block,ny_block,8), intent(out) :: & stPr ! stress combinations from replacement pressure @@ -1186,9 +1196,11 @@ subroutine calc_zeta_dPr (nx_block, ny_block, & capping = .false. - ! Initialize stPr and zetaD to zero (for cells where icetmask is false) - stPr = c0 - zetaD = c0 + ! Initialize stPr, zetax2 and etax2 to zero + ! (for cells where icetmask is false) + stPr = c0 + zetax2 = c0 + etax2 = c0 do ij = 1, icellt i = indxti(ij) @@ -1213,27 +1225,30 @@ subroutine calc_zeta_dPr (nx_block, ny_block, & Deltane , Deltanw , & Deltase , Deltasw) - if (capping) then - zetaD(i,j,1) = strength(i,j)/max(Deltane,tinyarea(i,j)) - zetaD(i,j,2) = strength(i,j)/max(Deltanw,tinyarea(i,j)) - zetaD(i,j,3) = strength(i,j)/max(Deltasw,tinyarea(i,j)) - zetaD(i,j,4) = strength(i,j)/max(Deltase,tinyarea(i,j)) - else - zetaD(i,j,1) = strength(i,j)/(Deltane + tinyarea(i,j)) - zetaD(i,j,2) = strength(i,j)/(Deltanw + tinyarea(i,j)) - zetaD(i,j,3) = strength(i,j)/(Deltasw + tinyarea(i,j)) - zetaD(i,j,4) = strength(i,j)/(Deltase + tinyarea(i,j)) - endif + !----------------------------------------------------------------- + ! viscous coefficients and replacement pressure + !----------------------------------------------------------------- + + call viscous_coeffs_and_rep_pressure (strength(i,j), tinyarea(i,j), & + Deltane, Deltanw, & + Deltasw, Deltase, & + zetax2(i,j,1), zetax2(i,j,2), & + zetax2(i,j,3), zetax2(i,j,4), & + etax2(i,j,1), etax2(i,j,2), & + etax2(i,j,3), etax2(i,j,4), & + rep_prs(i,j,1), rep_prs(i,j,2), & + rep_prs(i,j,3), rep_prs(i,j,4), & + capping) !----------------------------------------------------------------- ! the stresses ! kg/s^2 ! (1) northeast, (2) northwest, (3) southwest, (4) southeast !----------------------------------------------------------------- - stressp_1 = -zetaD(i,j,1)*(Deltane*(c1-Ktens)) - stressp_2 = -zetaD(i,j,2)*(Deltanw*(c1-Ktens)) - stressp_3 = -zetaD(i,j,3)*(Deltasw*(c1-Ktens)) - stressp_4 = -zetaD(i,j,4)*(Deltase*(c1-Ktens)) + stressp_1 = -rep_prs(i,j,1) + stressp_2 = -rep_prs(i,j,2) + stressp_3 = -rep_prs(i,j,3) + stressp_4 = -rep_prs(i,j,4) !----------------------------------------------------------------- ! combinations of the Pr related stresses for the momentum equation ! kg/s^2 @@ -1312,7 +1327,8 @@ subroutine stress_vp (nx_block , ny_block , & dxt , dyt , & cxp , cyp , & cxm , cym , & - zetaD , & + zetax2 , etax2 , & + rep_prs , & stressp_1 , stressp_2 , & stressp_3 , stressp_4 , & stressm_1 , stressm_2 , & @@ -1341,8 +1357,10 @@ subroutine stress_vp (nx_block , ny_block , & cxm ! 0.5*HTN - 1.5*HTS real (kind=dbl_kind), dimension(nx_block,ny_block,4), intent(in) :: & - zetaD ! 2*zeta - + zetax2 , & ! zetax2 = 2*zeta (bulk viscous coeff) + etax2 , & ! etax2 = 2*eta (shear viscous coeff) + rep_prs + real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: & stressp_1, stressp_2, stressp_3, stressp_4 , & ! sigma11+sigma22 stressm_1, stressm_2, stressm_3, stressm_4 , & ! sigma11-sigma22 @@ -1388,21 +1406,21 @@ subroutine stress_vp (nx_block , ny_block , & ! the stresses ! kg/s^2 ! (1) northeast, (2) northwest, (3) southwest, (4) southeast !----------------------------------------------------------------- - - stressp_1(i,j) = zetaD(i,j,1)*(divune*(c1+Ktens) - Deltane*(c1-Ktens)) - stressp_2(i,j) = zetaD(i,j,2)*(divunw*(c1+Ktens) - Deltanw*(c1-Ktens)) - stressp_3(i,j) = zetaD(i,j,3)*(divusw*(c1+Ktens) - Deltasw*(c1-Ktens)) - stressp_4(i,j) = zetaD(i,j,4)*(divuse*(c1+Ktens) - Deltase*(c1-Ktens)) - stressm_1(i,j) = zetaD(i,j,1)*tensionne*(c1+Ktens)*ecci - stressm_2(i,j) = zetaD(i,j,2)*tensionnw*(c1+Ktens)*ecci - stressm_3(i,j) = zetaD(i,j,3)*tensionsw*(c1+Ktens)*ecci - stressm_4(i,j) = zetaD(i,j,4)*tensionse*(c1+Ktens)*ecci + stressp_1(i,j) = zetax2(i,j,1)*divune - rep_prs(i,j,1) + stressp_2(i,j) = zetax2(i,j,2)*divunw - rep_prs(i,j,2) + stressp_3(i,j) = zetax2(i,j,3)*divusw - rep_prs(i,j,3) + stressp_4(i,j) = zetax2(i,j,4)*divuse - rep_prs(i,j,4) + + stressm_1(i,j) = etax2(i,j,1)*tensionne + stressm_2(i,j) = etax2(i,j,2)*tensionnw + stressm_3(i,j) = etax2(i,j,3)*tensionsw + stressm_4(i,j) = etax2(i,j,4)*tensionse - stress12_1(i,j) = zetaD(i,j,1)*shearne*p5*(c1+Ktens)*ecci - stress12_2(i,j) = zetaD(i,j,2)*shearnw*p5*(c1+Ktens)*ecci - stress12_3(i,j) = zetaD(i,j,3)*shearsw*p5*(c1+Ktens)*ecci - stress12_4(i,j) = zetaD(i,j,4)*shearse*p5*(c1+Ktens)*ecci + stress12_1(i,j) = etax2(i,j,1)*shearne*p5 + stress12_2(i,j) = etax2(i,j,2)*shearnw*p5 + stress12_3(i,j) = etax2(i,j,3)*shearsw*p5 + stress12_4(i,j) = etax2(i,j,4)*shearse*p5 enddo ! ij @@ -1534,7 +1552,7 @@ subroutine matvec (nx_block, ny_block, & cxm , cym , & uvel , vvel , & vrel , Cb , & - zetaD , & + zetax2 , etax2 , & umassdti, fm , & uarear , & Au , Av) @@ -1572,7 +1590,8 @@ subroutine matvec (nx_block, ny_block, & uarear ! 1/uarea real (kind=dbl_kind), dimension(nx_block,ny_block,4), intent(in) :: & - zetaD ! 2*zeta + zetax2 , & ! zetax2 = 2*zeta (bulk viscous coeff) + etax2 ! etax2 = 2*eta (shear viscous coeff) real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: & Au , & ! matvec, Fx = bx - Au (N/m^2) @@ -1647,20 +1666,20 @@ subroutine matvec (nx_block, ny_block, & ! NOTE: commented part of stressp is from the replacement 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)) + stressp_1 = zetax2(i,j,1)*divune! - Deltane*(c1-Ktens)) + stressp_2 = zetax2(i,j,2)*divunw! - Deltanw*(c1-Ktens)) + stressp_3 = zetax2(i,j,3)*divusw! - Deltasw*(c1-Ktens)) + stressp_4 = zetax2(i,j,4)*divuse! - 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 + stressm_1 = etax2(i,j,1)*tensionne + stressm_2 = etax2(i,j,2)*tensionnw + stressm_3 = etax2(i,j,3)*tensionsw + stressm_4 = etax2(i,j,4)*tensionse - 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 + stress12_1 = etax2(i,j,1)*shearne*p5 + stress12_2 = etax2(i,j,2)*shearnw*p5 + stress12_3 = etax2(i,j,3)*shearsw*p5 + stress12_4 = etax2(i,j,4)*shearse*p5 !----------------------------------------------------------------- ! combinations of the stresses for the momentum equation ! kg/s^2 @@ -1991,7 +2010,8 @@ subroutine formDiag_step1 (nx_block, ny_block, & dxhy , dyhx , & cxp , cyp , & cxm , cym , & - zetaD , Drheo) + zetax2 , etax2 , & + Drheo) integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions @@ -2012,7 +2032,8 @@ subroutine formDiag_step1 (nx_block, ny_block, & cxm ! 0.5*HTN - 1.5*HTS real (kind=dbl_kind), dimension(nx_block,ny_block,4), intent(in) :: & - zetaD ! 2*zeta + zetax2 , & ! zetax2 = 2*zeta (bulk viscous coeff) + etax2 ! etax2 = 2*eta (shear viscous coeff) real (kind=dbl_kind), dimension(nx_block,ny_block,8), intent(out) :: & Drheo ! intermediate value for diagonal components of matrix A associated @@ -2200,20 +2221,20 @@ subroutine formDiag_step1 (nx_block, ny_block, & ! (1) northeast, (2) northwest, (3) southwest, (4) southeast !----------------------------------------------------------------- - stressp_1 = zetaD(i,j,1)*divune*(c1+Ktens) - stressp_2 = zetaD(i,j,2)*divunw*(c1+Ktens) - stressp_3 = zetaD(i,j,3)*divusw*(c1+Ktens) - stressp_4 = zetaD(i,j,4)*divuse*(c1+Ktens) + stressp_1 = zetax2(i,j,1)*divune + stressp_2 = zetax2(i,j,2)*divunw + stressp_3 = zetax2(i,j,3)*divusw + stressp_4 = zetax2(i,j,4)*divuse - 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 + stressm_1 = etax2(i,j,1)*tensionne + stressm_2 = etax2(i,j,2)*tensionnw + stressm_3 = etax2(i,j,3)*tensionsw + stressm_4 = etax2(i,j,4)*tensionse - 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 + stress12_1 = etax2(i,j,1)*shearne*p5 + stress12_2 = etax2(i,j,2)*shearnw*p5 + stress12_3 = etax2(i,j,3)*shearsw*p5 + stress12_4 = etax2(i,j,4)*shearse*p5 !----------------------------------------------------------------- ! combinations of the stresses for the momentum equation ! kg/s^2 @@ -2657,7 +2678,7 @@ end subroutine qr_delete ! ! authors: Stéphane Gaudreault, Abdessamad Qaddouri, Philippe Blain, ECCC - subroutine fgmres (zetaD , & + subroutine fgmres (zetax2 , etax2 , & Cb , vrel , & umassdti , & halo_info_mask , & @@ -2673,8 +2694,9 @@ subroutine fgmres (zetaD , & use ice_timers, only: ice_timer_start, ice_timer_stop, timer_bound real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks,4), intent(in) :: & - zetaD ! zetaD = 2*zeta (viscous coefficient) - + zetax2 , & ! zetax2 = 2*zeta (bulk viscous coeff) + etax2 ! etax2 = 2*eta (shear viscous coeff) + real (kind=dbl_kind), dimension(nx_block, ny_block, max_blocks), intent(in) :: & vrel , & ! coefficient for tauw Cb , & ! seabed stress coefficient @@ -2784,7 +2806,7 @@ subroutine fgmres (zetaD , & cxm (:,:,iblk) , cym (:,:,iblk), & solx (:,:,iblk) , soly (:,:,iblk), & vrel (:,:,iblk) , Cb (:,:,iblk), & - zetaD (:,:,iblk,:), & + zetax2 (:,:,iblk,:), etax2 (:,:,iblk,:), & umassdti (:,:,iblk) , fm (:,:,iblk), & uarear (:,:,iblk) , & workspace_x(:,:,iblk) , workspace_y(:,:,iblk)) @@ -2855,7 +2877,7 @@ subroutine fgmres (zetaD , & initer = initer + 1 nextit = initer + 1 ! precondition the current Arnoldi vector - call precondition(zetaD , & + call precondition(zetax2 , etax2 , & Cb , vrel , & umassdti , & arnoldi_basis_x(:,:,:,initer), & @@ -2891,7 +2913,7 @@ subroutine fgmres (zetaD , & cxm (:,:,iblk) , cym (:,:,iblk), & workspace_x(:,:,iblk) , workspace_y(:,:,iblk), & vrel (:,:,iblk) , Cb (:,:,iblk), & - zetaD (:,:,iblk,:), & + zetax2 (:,:,iblk,:), etax2 (:,:,iblk,:), & umassdti (:,:,iblk) , fm (:,:,iblk), & uarear (:,:,iblk) , & arnoldi_basis_x(:,:,iblk,nextit), & @@ -3057,7 +3079,7 @@ end subroutine fgmres ! ! authors: Stéphane Gaudreault, Abdessamad Qaddouri, Philippe Blain, ECCC - subroutine pgmres (zetaD , & + subroutine pgmres (zetax2 , etax2 , & Cb , vrel , & umassdti , & bx , by , & @@ -3068,8 +3090,9 @@ subroutine pgmres (zetaD , & nbiter) real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks,4), intent(in) :: & - zetaD ! zetaD = 2*zeta (viscous coefficient) - + zetax2 , & ! zetax2 = 2*zeta (bulk viscous coeff) + etax2 ! etax2 = 2*eta (shear viscous coeff) + real (kind=dbl_kind), dimension(nx_block, ny_block, max_blocks), intent(in) :: & vrel , & ! coefficient for tauw Cb , & ! seabed stress coefficient @@ -3176,7 +3199,7 @@ subroutine pgmres (zetaD , & cxm (:,:,iblk) , cym (:,:,iblk), & solx (:,:,iblk) , soly (:,:,iblk), & vrel (:,:,iblk) , Cb (:,:,iblk), & - zetaD (:,:,iblk,:), & + zetax2 (:,:,iblk,:), etax2 (:,:,iblk,:), & umassdti (:,:,iblk) , fm (:,:,iblk), & uarear (:,:,iblk) , & workspace_x(:,:,iblk) , workspace_y(:,:,iblk)) @@ -3248,7 +3271,7 @@ subroutine pgmres (zetaD , & nextit = initer + 1 ! precondition the current Arnoldi vector - call precondition(zetaD , & + call precondition(zetax2 , etax2 , & Cb , vrel , & umassdti , & arnoldi_basis_x(:,:,:,initer), & @@ -3272,7 +3295,7 @@ subroutine pgmres (zetaD , & cxm (:,:,iblk) , cym (:,:,iblk), & workspace_x(:,:,iblk) , workspace_y(:,:,iblk), & vrel (:,:,iblk) , Cb (:,:,iblk), & - zetaD (:,:,iblk,:), & + zetax2 (:,:,iblk,:), etax2 (:,:,iblk,:), & umassdti (:,:,iblk) , fm (:,:,iblk), & uarear (:,:,iblk) , & arnoldi_basis_x(:,:,iblk,nextit), & @@ -3385,7 +3408,7 @@ subroutine pgmres (zetaD , & end do ! Call preconditioner - call precondition(zetaD , & + call precondition(zetax2 , etax2 , & Cb , vrel , & umassdti , & workspace_x , workspace_y, & @@ -3451,7 +3474,7 @@ end subroutine pgmres ! ! authors: Philippe Blain, ECCC - subroutine precondition(zetaD , & + subroutine precondition(zetax2 , etax2, & Cb , vrel , & umassdti , & vx , vy , & @@ -3460,8 +3483,9 @@ subroutine precondition(zetaD , & wx , wy) real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks,4), intent(in) :: & - zetaD ! zetaD = 2*zeta (viscous coefficient) - + zetax2 , & ! zetax2 = 2*zeta (bulk viscous coeff) + etax2 ! etax2 = 2*eta (shear viscous coeff) + real (kind=dbl_kind), dimension(nx_block, ny_block, max_blocks), intent(in) :: & vrel , & ! coefficient for tauw Cb , & ! seabed stress coefficient @@ -3525,7 +3549,7 @@ subroutine precondition(zetaD , & tolerance = reltol_pgmres maxinner = dim_pgmres maxouter = maxits_pgmres - call pgmres (zetaD, & + call pgmres (zetax2, etax2 , & Cb , vrel , & umassdti , & vx , vy , &