diff --git a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 index 42fe1b53d..b58fbc1c7 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90 @@ -590,8 +590,7 @@ subroutine evp (dt) if (seabed_stress) then - select case (trim(grid_ice)) - case('B') + if (grid_ice == "B") then if ( seabed_stress_method == 'LKD' ) then !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime) @@ -616,7 +615,7 @@ subroutine evp (dt) !$OMP END PARALLEL DO endif - case('CD','C') + elseif (grid_ice == "C" .or. grid_ice == "CD") then if ( seabed_stress_method == 'LKD' ) then !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime) @@ -649,7 +648,7 @@ subroutine evp (dt) !$OMP END PARALLEL DO endif - end select + endif endif @@ -691,8 +690,7 @@ subroutine evp (dt) do ksub = 1,ndte ! subcycling - select case (grid_ice) - case('B') + if (grid_ice == "B") then !$OMP PARALLEL DO PRIVATE(iblk,strtmp) SCHEDULE(runtime) do iblk = 1, nblocks @@ -757,7 +755,7 @@ subroutine evp (dt) enddo ! iblk !$OMP END PARALLEL DO - case('CD','C') + elseif (grid_ice == "C" .or. grid_ice == "CD") then !$OMP PARALLEL DO PRIVATE(iblk) do iblk = 1, nblocks @@ -967,7 +965,7 @@ subroutine evp (dt) uvel(:,:,:) = uvel(:,:,:)*uvm(:,:,:) vvel(:,:,:) = vvel(:,:,:)*uvm(:,:,:) - end select + endif ! grid_ice call ice_timer_start(timer_bound) call stack_velocity_field(uvel, vvel, fld2) @@ -1169,9 +1167,8 @@ subroutine stress (nx_block, ny_block, & stress12_3, stress12_4, & str ) - use ice_dyn_shared, only: strain_rates, viscous_coeffs_and_rep_pressure_T, & - capping - + use ice_dyn_shared, only: strain_rates, visccoeff_replpress, capping + integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions icellt ! no. of cells where icetmask = 1 @@ -1262,27 +1259,18 @@ subroutine stress (nx_block, ny_block, & ! viscous coefficients and replacement pressure !----------------------------------------------------------------- - call viscous_coeffs_and_rep_pressure_T (strength(i,j), DminTarea(i,j),& - Deltane, zetax2ne, & - etax2ne, rep_prsne, & - capping) - - call viscous_coeffs_and_rep_pressure_T (strength(i,j), DminTarea(i,j),& - Deltanw, zetax2nw, & - etax2nw, rep_prsnw, & - capping) + call visccoeff_replpress (strength(i,j), DminTarea(i,j), Deltane, & + zetax2ne, etax2ne, rep_prsne, capping) - call viscous_coeffs_and_rep_pressure_T (strength(i,j), DminTarea(i,j),& - Deltasw, zetax2sw, & - etax2sw, rep_prssw, & - capping) + call visccoeff_replpress (strength(i,j), DminTarea(i,j), Deltanw, & + zetax2nw, etax2nw, rep_prsnw, capping) - call viscous_coeffs_and_rep_pressure_T (strength(i,j), DminTarea(i,j),& - Deltase, zetax2se, & - etax2se, rep_prsse, & - capping) + call visccoeff_replpress (strength(i,j), DminTarea(i,j), Deltasw, & + zetax2sw, etax2sw, rep_prssw, capping) + + call visccoeff_replpress (strength(i,j), DminTarea(i,j), Deltase, & + zetax2se, etax2se, rep_prsse, capping) - !----------------------------------------------------------------- ! the stresses ! kg/s^2 ! (1) northeast, (2) northwest, (3) southwest, (4) southeast @@ -1471,7 +1459,7 @@ subroutine stress_T (nx_block, ny_block, & stress12T ) use ice_dyn_shared, only: strain_rates_T, capping, & - viscous_coeffs_and_rep_pressure_T + visccoeff_replpress integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions @@ -1538,10 +1526,10 @@ subroutine stress_T (nx_block, ny_block, & ! viscous coefficients and replacement pressure at T point !----------------------------------------------------------------- - call viscous_coeffs_and_rep_pressure_T (strength(i,j), & - DminTarea(i,j), DeltaT, & - zetax2T(i,j),etax2T(i,j),& - rep_prsT, capping ) + call visccoeff_replpress (strength(i,j), DminTarea(i,j), & + DeltaT , zetax2T (i,j), & + etax2T (i,j), rep_prsT , & + capping) !----------------------------------------------------------------- ! the stresses ! kg/s^2 @@ -1587,10 +1575,10 @@ subroutine stress_U (nx_block, ny_block, & stress12U ) use ice_dyn_shared, only: strain_rates_U, & - viscous_coeffs_and_rep_pressure_T2U, & - viscous_coeffs_and_rep_pressure_U, & + visccoeff_replpress_avgstr, & + visccoeff_replpress_avgzeta, & visc_coeff_method, deltaminEVP, capping - + integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions icellu ! no. of cells where iceumask = 1 @@ -1637,7 +1625,7 @@ subroutine stress_U (nx_block, ny_block, & real (kind=dbl_kind) :: & divU, tensionU, shearU, DeltaU, & ! strain rates at U point zetax2U, etax2U, rep_prsU, & ! replacement pressure at U point - DminUarea + DminUarea, strtmp, areatmp ! Dmin on U and tmp variables character(len=*), parameter :: subname = '(stress_U)' @@ -1668,28 +1656,26 @@ subroutine stress_U (nx_block, ny_block, & !----------------------------------------------------------------- if (visc_coeff_method == 'avg_zeta') then - call viscous_coeffs_and_rep_pressure_T2U (zetax2T(i ,j ), zetax2T(i ,j+1), & - zetax2T(i+1,j+1), zetax2T(i+1,j ), & - etax2T (i ,j ), etax2T (i ,j+1), & - etax2T (i+1,j+1), etax2T (i+1,j ), & - hm (i ,j ), hm (i ,j+1), & - hm (i+1,j+1), hm (i+1,j ), & - tarea (i ,j ), tarea (i ,j+1), & - tarea (i+1,j+1), tarea (i+1,j ), & - DeltaU,zetax2U, etax2U, rep_prsU) - elseif (visc_coeff_method == 'avg_strength') then + call visccoeff_replpress_avgzeta (zetax2T (i ,j ), zetax2T (i ,j+1), & + zetax2T (i+1,j+1), zetax2T (i+1,j ), & + etax2T (i ,j ), etax2T (i ,j+1), & + etax2T (i+1,j+1), etax2T (i+1,j ), & + hm (i ,j ), hm (i ,j+1), & + hm (i+1,j+1), hm (i+1,j ), & + tarea (i ,j ), tarea (i ,j+1), & + tarea (i+1,j+1), tarea (i+1,j ), & + DeltaU, zetax2U, etax2U, rep_prsU) + elseif (visc_coeff_method == 'avg_strength') then DminUarea = deltaminEVP*uarea(i,j) - call viscous_coeffs_and_rep_pressure_U (strength(i ,j ), strength(i ,j+1), & - strength(i+1,j+1), strength(i+1,j ), & - hm (i ,j ) , hm (i ,j+1), & - hm (i+1,j+1) , hm (i+1,j ), & - tarea (i ,j ) , tarea (i ,j+1), & - tarea (i+1,j+1) , tarea (i+1,j ), & - DminUarea, & - DeltaU , capping, & - zetax2U, etax2U, rep_prsU) - + call visccoeff_replpress_avgstr (strength(i ,j ), strength(i ,j+1), & + strength(i+1,j+1), strength(i+1,j ), & + hm (i ,j ) , hm (i ,j+1), & + hm (i+1,j+1) , hm (i+1,j ), & + tarea (i ,j ) , tarea (i ,j+1), & + tarea (i+1,j+1) , tarea (i+1,j ), & + DminUarea, DeltaU, & + zetax2U, etax2U, rep_prsU, capping) endif !----------------------------------------------------------------- @@ -1779,8 +1765,7 @@ subroutine div_stress (nx_block, ny_block, & ! F1,F2 : div of stress tensor for u,v components !----------------------------------------------------------------- - select case (trim(grid_location)) - case('E') + if (grid_location == "E") then F1(i,j) = arear(i,j) * & ( p5 * dyE_N(i,j) * ( stresspF1(i+1,j)-stresspF1(i,j) ) & @@ -1796,7 +1781,7 @@ subroutine div_stress (nx_block, ny_block, & + (c1/dyE_N(i,j)) * ( (dyT_U(i+1,j)**2) * stress12F2(i+1,j) & -(dyT_U(i,j)**2)*stress12F2(i,j) ) ) - case('N') + elseif (grid_location == "N") then F1(i,j) = arear(i,j) * & ( p5 * dyE_N(i,j) * ( stresspF1(i,j)-stresspF1(i-1,j) ) & @@ -1811,9 +1796,9 @@ subroutine div_stress (nx_block, ny_block, & -(dxT_U(i,j)**2)*stressmF2(i,j) ) & + (c1/dyE_N(i,j)) * ( (dyT_U(i,j)**2) * stress12F2(i,j) & -(dyT_U(i-1,j)**2)*stress12F2(i-1,j) ) ) - case default + else call abort_ice(subname // ' unknown grid_location: ' // grid_location) - end select + endif enddo ! ij diff --git a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 index 8301a967b..f9f1b3c92 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_shared.F90 @@ -29,10 +29,9 @@ module ice_dyn_shared alloc_dyn_shared, & deformations, deformations_T, & strain_rates, strain_rates_T, strain_rates_U, & - viscous_coeffs_and_rep_pressure, & - viscous_coeffs_and_rep_pressure_T, & - viscous_coeffs_and_rep_pressure_T2U, & - viscous_coeffs_and_rep_pressure_U, & + visccoeff_replpress, & + visccoeff_replpress_avgstr, & + visccoeff_replpress_avgzeta, & stack_velocity_field, unstack_velocity_field ! namelist parameters @@ -52,6 +51,10 @@ module ice_dyn_shared character (len=char_len), public :: & evp_algorithm ! standard_2d = 2D org version (standard) ! shared_mem_1d = 1d without mpi call and refactorization to 1d + + real (kind=dbl_kind), public :: & + elasticDamp ! coefficient for calculating the parameter E, elastic damping parameter + ! other EVP parameters character (len=char_len), public :: & @@ -62,7 +65,6 @@ module ice_dyn_shared real (kind=dbl_kind), parameter, public :: & - eyc = 0.36_dbl_kind, & ! coefficient for calculating the parameter E u0 = 5e-5_dbl_kind, & ! residual velocity for seabed stress (m/s) cosw = c1 , & ! cos(ocean turning angle) ! turning angle = 0 sinw = c0 , & ! sin(ocean turning angle) ! turning angle = 0 @@ -190,7 +192,7 @@ subroutine init_dyn (dt) if (my_task == master_task) then write(nu_diag,*) 'dt = ',dt write(nu_diag,*) 'dte = ',dt/real(ndte,kind=dbl_kind) - write(nu_diag,*) 'tdamp =', eyc*dt + write(nu_diag,*) 'tdamp =', elasticDamp * dt endif allocate(fcor_blk(nx_block,ny_block,max_blocks)) @@ -318,8 +320,8 @@ subroutine set_evp_parameters (dt) ecci = c1/e_yieldcurve**2 ! temporary for 1d evp ! constants for stress equation - !tdamp2 = c2*eyc*dt ! s - !dte2T = dte/tdamp2 or c1/(c2*eyc*real(ndte,kind=dbl_kind)) ! ellipse (unitless) + !tdamp2 = c2 * elasticDamp * dt ! s + !dte2T = dte/tdamp2 or c1/(c2*elasticDamp*real(ndte,kind=dbl_kind)) ! ellipse (unitless) if (revised_evp) then ! Bouillon et al, Ocean Mod 2013 revp = c1 @@ -330,7 +332,7 @@ subroutine set_evp_parameters (dt) !arlx1i = dte2T !arlx = c1/arlx1i !brlx = dt*dtei - arlx = c2*eyc*real(ndte,kind=dbl_kind) + arlx = c2 * elasticDamp * real(ndte,kind=dbl_kind) arlx1i = c1/arlx brlx = real(ndte,kind=dbl_kind) denom1 = c1/(c1+arlx1i) @@ -1032,7 +1034,7 @@ subroutine stepu_Cgrid (nx_block, ny_block, & enddo ! ij - end subroutine stepu_Cgrid + end subroutine stepu_Cgrid !======================================================================= @@ -1132,7 +1134,7 @@ subroutine stepv_Cgrid (nx_block, ny_block, & enddo ! ij - end subroutine stepv_Cgrid + end subroutine stepv_Cgrid !======================================================================= @@ -1306,7 +1308,7 @@ subroutine seabed_stress_factor_LKD (nx_block, ny_block, & enddo ! ij - end subroutine seabed_stress_factor_LKD + end subroutine seabed_stress_factor_LKD !======================================================================= ! Computes seabed (basal) stress factor Tbu (landfast ice) based on @@ -1491,39 +1493,38 @@ subroutine seabed_stress_factor_prob (nx_block, ny_block, & endif enddo - select case (trim(grid_ice)) - case('B') - do ij = 1, icellu - i = indxui(ij) - j = indxuj(ij) - ! convert quantities to U-location - Tbu(i,j) = grid_neighbor_max(Tbt, i, j, 'U') - enddo ! ij - case('CD','C') - if(present(Tbe) .and. present(TbN) .and. & - present(icelle) .and. present(icelln) .and. & - present(indxei) .and. present(indxej) .and. & - present(indxni) .and. present(indxnj)) then - - do ij = 1, icelle - i = indxei(ij) - j = indxej(ij) - ! convert quantities to E-location - TbE(i,j) = grid_neighbor_max(Tbt, i, j, 'E') - enddo - do ij = 1, icelln - i = indxni(ij) - j = indxnj(ij) - ! convert quantities to N-location - TbN(i,j) = grid_neighbor_max(Tbt, i, j, 'N') - enddo + if (grid_ice == "B") then + do ij = 1, icellu + i = indxui(ij) + j = indxuj(ij) + ! convert quantities to U-location + Tbu(i,j) = grid_neighbor_max(Tbt, i, j, 'U') + enddo ! ij + elseif (grid_ice == "C" .or. grid_ice == "CD") then + if (present(Tbe) .and. present(TbN) .and. & + present(icelle) .and. present(icelln) .and. & + present(indxei) .and. present(indxej) .and. & + present(indxni) .and. present(indxnj)) then + + do ij = 1, icelle + i = indxei(ij) + j = indxej(ij) + ! convert quantities to E-location + TbE(i,j) = grid_neighbor_max(Tbt, i, j, 'E') + enddo + do ij = 1, icelln + i = indxni(ij) + j = indxnj(ij) + ! convert quantities to N-location + TbN(i,j) = grid_neighbor_max(Tbt, i, j, 'N') + enddo - else - call abort_ice(subname // ' insufficient number of arguments for grid_ice:' // grid_ice) - endif - end select + else + call abort_ice(subname // ' insufficient number of arguments for grid_ice:' // grid_ice) + endif + endif - end subroutine seabed_stress_factor_prob + end subroutine seabed_stress_factor_prob !======================================================================= @@ -1688,7 +1689,7 @@ subroutine deformations (nx_block, ny_block, & enddo ! ij - end subroutine deformations + end subroutine deformations !======================================================================= @@ -1780,7 +1781,7 @@ subroutine deformations_T (nx_block, ny_block, & enddo ! ij - end subroutine deformations_T + end subroutine deformations_T !======================================================================= @@ -1929,9 +1930,8 @@ subroutine strain_rates_T (nx_block, ny_block, & ! Delta (in the denominator of zeta, eta) DeltaT = sqrt(divT**2 + e_factor*(tensionT**2 + shearT**2)) - end subroutine strain_rates_T + end subroutine strain_rates_T - !======================================================================= ! Compute strain rates at the U point including boundary conditions @@ -2034,100 +2034,27 @@ subroutine strain_rates_U (nx_block, ny_block, & ! Delta (in the denominator of zeta, eta) DeltaU = sqrt(divU**2 + e_factor*(tensionU**2 + shearU**2)) - end subroutine strain_rates_U - - !======================================================================= - ! Computes viscous coefficients and replacement pressure for stress - ! calculations. Note that tensile strength is included here. - ! - ! Hibler, W. D. (1979). A dynamic thermodynamic sea ice model. J. Phys. - ! Oceanogr., 9, 817-846. - ! - ! Konig Beatty, C. and Holland, D. M. (2010). Modeling landfast ice by - ! adding tensile strength. J. Phys. Oceanogr. 40, 185-198. - ! - ! Lemieux, J. F. et al. (2016). Improving the simulation of landfast ice - ! by combining tensile strength and a parameterization for grounded ridges. - ! J. Geophys. Res. Oceans, 121, 7354-7368. - - subroutine viscous_coeffs_and_rep_pressure (strength, DminTarea,& - Deltane, Deltanw, & - Deltasw, Deltase, & - zetax2ne, zetax2nw, & - zetax2sw, zetax2se, & - etax2ne, etax2nw, & - etax2sw, etax2se, & - rep_prsne, rep_prsnw,& - rep_prssw, rep_prsse,& - capping) - - real (kind=dbl_kind), intent(in):: & - strength, DminTarea ! at the t-point - - real (kind=dbl_kind), intent(in):: & - Deltane, Deltanw, Deltasw, Deltase ! Delta at each corner - - real(kind=dbl_kind) , intent(in):: capping - - real (kind=dbl_kind), intent(out):: & - 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) :: & - tmpcalcne, tmpcalcnw, tmpcalcsw, tmpcalcse + end subroutine strain_rates_U - ! NOTE: for comp. efficiency 2 x zeta and 2 x eta are used in the code - - tmpcalcne = capping *(strength/max(Deltane, DminTarea))+ & - (c1-capping)* strength/ (Deltane+ DminTarea) - tmpcalcnw = capping *(strength/max(Deltanw, DminTarea))+ & - (c1-capping)* strength/ (Deltanw+ DminTarea) - tmpcalcsw = capping *(strength/max(Deltasw, DminTarea))+ & - (c1-capping)* strength/ (Deltasw+ DminTarea) - tmpcalcse = capping *(strength/max(Deltase, DminTarea))+ & - (c1-capping)* strength/ (Deltase+ DminTarea) - - zetax2ne = (c1+Ktens)*tmpcalcne ! northeast - rep_prsne = (c1-Ktens)*tmpcalcne*Deltane - etax2ne = epp2i*zetax2ne - - zetax2nw = (c1+Ktens)*tmpcalcnw ! northwest - rep_prsnw = (c1-Ktens)*tmpcalcnw*Deltanw - etax2nw = epp2i*zetax2nw +!======================================================================= +! Computes viscous coefficients and replacement pressure for stress +! calculations. Note that tensile strength is included here. +! +! Hibler, W. D. (1979). A dynamic thermodynamic sea ice model. J. Phys. +! Oceanogr., 9, 817-846. +! +! Konig Beatty, C. and Holland, D. M. (2010). Modeling landfast ice by +! adding tensile strength. J. Phys. Oceanogr. 40, 185-198. +! +! Lemieux, J. F. et al. (2016). Improving the simulation of landfast ice +! by combining tensile strength and a parameterization for grounded ridges. +! J. Geophys. Res. Oceans, 121, 7354-7368. - zetax2sw = (c1+Ktens)*tmpcalcsw ! southwest - rep_prssw = (c1-Ktens)*tmpcalcsw*Deltasw - etax2sw = epp2i*zetax2sw - - zetax2se = (c1+Ktens)*tmpcalcse ! southeast - rep_prsse = (c1-Ktens)*tmpcalcse*Deltase - etax2se = epp2i*zetax2se - - end subroutine viscous_coeffs_and_rep_pressure - - !======================================================================= - ! Computes viscous coefficients and replacement pressure for stress - ! calculations. Note that tensile strength is included here. - ! - ! Hibler, W. D. (1979). A dynamic thermodynamic sea ice model. J. Phys. - ! Oceanogr., 9, 817-846. - ! - ! Konig Beatty, C. and Holland, D. M. (2010). Modeling landfast ice by - ! adding tensile strength. J. Phys. Oceanogr. 40, 185-198. - ! - ! Lemieux, J. F. et al. (2016). Improving the simulation of landfast ice - ! by combining tensile strength and a parameterization for grounded ridges. - ! J. Geophys. Res. Oceans, 121, 7354-7368. - - subroutine viscous_coeffs_and_rep_pressure_T (strength, DminTarea, & - Delta , zetax2 , & - etax2 , rep_prs , & - capping) + subroutine visccoeff_replpress(strength, DminArea, Delta, & + zetax2, etax2, rep_prs, capping) real (kind=dbl_kind), intent(in):: & - strength, DminTarea + strength, DminArea real (kind=dbl_kind), intent(in):: & Delta, capping @@ -2139,123 +2066,115 @@ subroutine viscous_coeffs_and_rep_pressure_T (strength, DminTarea, & real (kind=dbl_kind) :: & tmpcalc - character(len=*), parameter :: subname = '(viscous_coeffs_and_rep_pressure_T)' + character(len=*), parameter :: subname = '(visccoeff_replpress)' ! NOTE: for comp. efficiency 2 x zeta and 2 x eta are used in the code - tmpcalc = capping *(strength/max(Delta,DminTarea))+ & - (c1-capping)*(strength/(Delta + DminTarea)) - zetax2 = (c1+Ktens)*tmpcalc + tmpcalc = capping *(strength/max(Delta,DminArea))+ & + (c1-capping)*(strength/(Delta + DminArea)) + zetax2 = (c1+Ktens)*tmpcalc rep_prs = (c1-Ktens)*tmpcalc*Delta - etax2 = epp2i*zetax2 + etax2 = epp2i*zetax2 - end subroutine viscous_coeffs_and_rep_pressure_T + end subroutine visccoeff_replpress +!======================================================================= - subroutine viscous_coeffs_and_rep_pressure_T2U (zetax2T_00, zetax2T_01, & - zetax2T_11, zetax2T_10, & - etax2T_00, etax2T_01, & - etax2T_11, etax2T_10, & - maskT_00, maskT_01, & - maskT_11, maskT_10, & - tarea_00, tarea_01, & - tarea_11, tarea_10, & - deltaU, & - zetax2U, etax2U, & - rep_prsU) - + subroutine visccoeff_replpress_avgzeta (zetax2T1, zetax2T2, & + zetax2T3, zetax2T4, & + etax2T1, etax2T2, & + etax2T3, etax2T4, & + mask1, mask2, & + mask3, mask4, & + area1, area2, & + area3, area4, & + deltaU, zetax2U, etax2U, rep_prsU) real (kind=dbl_kind), intent(in):: & - zetax2T_00,zetax2T_10,zetax2T_11,zetax2T_01, & - etax2T_00, etax2T_10, etax2T_11, etax2T_01, & ! 2 x viscous coeffs, replacement pressure - maskT_00, maskT_10, maskT_11, maskT_01, & - tarea_00, tarea_10, tarea_11, tarea_01, & + zetax2T1,zetax2T2,zetax2T3,zetax2T4, & + etax2T1, etax2T2, etax2T3, etax2T4, & ! 2 x viscous coeffs, replacement pressure + mask1, mask2, mask3, mask4, & + area1, area2, area3, area4, & deltaU real (kind=dbl_kind), intent(out):: zetax2U, etax2U, rep_prsU ! local variables - + real (kind=dbl_kind) :: & - Totarea + areatmp - character(len=*), parameter :: subname = '(viscous_coeffs_and_rep_pressure_T2U)' + character(len=*), parameter :: subname = '(visccoeff_replpress_avgzeta)' ! NOTE: for comp. efficiency 2 x zeta and 2 x eta are used in the code - Totarea = maskT_00*Tarea_00 + & - maskT_10*Tarea_10 + & - maskT_11*Tarea_11 + & - maskT_01*Tarea_01 - zetax2U = (maskT_00*Tarea_00 *zetax2T_00 + & - maskT_10*Tarea_10 *zetax2T_10 + & - maskT_11*Tarea_11 *zetax2T_11 + & - maskT_01*Tarea_01 *zetax2T_01)/Totarea - - etax2U = (maskT_00*Tarea_00 *etax2T_00 + & - maskT_10*Tarea_10 *etax2T_10 + & - maskT_11*Tarea_11 *etax2T_11 + & - maskT_01*Tarea_01 *etax2T_01)/Totarea - + + areatmp = (mask1 * area1 + & + mask4 * area4 + & + mask3 * area3 + & + mask2 * area2) + + zetax2U = (mask1 * area1 * zetax2T1 + & + mask4 * area4 * zetax2T4 + & + mask3 * area3 * zetax2T3 + & + mask2 * area2 * zetax2T2) / areatmp + + etax2U = (mask1 * area1 * etax2T1 + & + mask4 * area4 * etax2T4 + & + mask3 * area3 * etax2T3 + & + mask2 * area2 * etax2T2) / areatmp + rep_prsU = (c1-Ktens)/(c1+Ktens)*zetax2U*deltaU - end subroutine viscous_coeffs_and_rep_pressure_T2U + end subroutine visccoeff_replpress_avgzeta !======================================================================= - subroutine viscous_coeffs_and_rep_pressure_U (strength_00, strength_01, & - strength_11, strength_10, & - maskT_00, maskT_01, & - maskT_11, maskT_10, & - tarea_00, tarea_01, & - tarea_11, tarea_10, & - DminUarea, & - deltaU, capping, & - zetax2U, etax2U, & - rep_prsU) - + subroutine visccoeff_replpress_avgstr (strength1, strength2, & + strength3, strength4, & + mask1, mask2, & + mask3, mask4, & + area1, area2, & + area3, area4, & + DminUarea, deltaU, & + zetax2U, etax2U, rep_prsU, capping) real (kind=dbl_kind), intent(in):: & - strength_00,strength_10,strength_11,strength_01, & - maskT_00, maskT_10, maskT_11, maskT_01, & - tarea_00, tarea_10, tarea_11, tarea_01, & + strength1,strength2,strength3,strength4, & + mask1, mask2, mask3, mask4, & + area1, area2, area3, area4, & DminUarea, deltaU, capping real (kind=dbl_kind), intent(out):: zetax2U, etax2U, rep_prsU ! local variables - + real (kind=dbl_kind) :: & - Totarea, tmpcalc, strength + areatmp, strtmp ! area and strength average - character(len=*), parameter :: subname = '(viscous_coeffs_and_rep_pressure_U)' + character(len=*), parameter :: subname = '(visccoeff_replpress_avgstr)' ! NOTE: for comp. efficiency 2 x zeta and 2 x eta are used in the code - Totarea = maskT_00*Tarea_00 + & - maskT_10*Tarea_10 + & - maskT_11*Tarea_11 + & - maskT_01*Tarea_01 - strength = (maskT_00*Tarea_00 *strength_00 + & - maskT_10*Tarea_10 *strength_10 + & - maskT_11*Tarea_11 *strength_11 + & - maskT_01*Tarea_01 *strength_01)/Totarea - - ! rep_prsU = (c1-Ktens)/(c1+Ktens)*zetax2U*deltaU - ! IMPROVE the calc below are the same as in the other viscous coeff...could reduce redundency - ! we could have a strength_U subroutine and then calc the visc coeff - - tmpcalc = capping *(strength/max(deltaU,DminUarea))+ & - (c1-capping)*(strength/(deltaU + DminUarea)) - zetax2U = (c1+Ktens)*tmpcalc - rep_prsU = (c1-Ktens)*tmpcalc*deltaU - etax2U = epp2i*zetax2U + + areatmp = (mask1 * area1 + & + mask4 * area4 + & + mask3 * area3 + & + mask2 * area2) + + strtmp = (mask1 * area1 * strength1 + & + mask4 * area4 * strength4 + & + mask3 * area3 * strength3 + & + mask2 * area2 * strength2) / areatmp + + call visccoeff_replpress (strtmp, DminUarea, deltaU, & + zetax2U, etax2U, rep_prsU, capping) - end subroutine viscous_coeffs_and_rep_pressure_U + end subroutine visccoeff_replpress_avgstr !======================================================================= ! Load velocity components into array for boundary updates - subroutine stack_velocity_field(uvel, vvel, fld2) + subroutine stack_velocity_field(uvel, vvel, fld2) use ice_domain, only: nblocks @@ -2287,7 +2206,7 @@ end subroutine stack_velocity_field ! Unload velocity components from array after boundary updates - subroutine unstack_velocity_field(fld2, uvel, vvel) + subroutine unstack_velocity_field(fld2, uvel, vvel) use ice_domain, only: nblocks diff --git a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 index 5d414c204..17de3bca8 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 @@ -1157,8 +1157,8 @@ subroutine calc_zeta_dPr (nx_block, ny_block, & zetax2 , etax2 , & rep_prs , stPr) - use ice_dyn_shared, only: strain_rates, viscous_coeffs_and_rep_pressure, & - capping + use ice_dyn_shared, only: strain_rates, visccoeff_replpress, & + capping integer (kind=int_kind), intent(in) :: & nx_block, ny_block, & ! block dimensions @@ -1240,17 +1240,23 @@ subroutine calc_zeta_dPr (nx_block, ny_block, & ! viscous coefficients and replacement pressure !----------------------------------------------------------------- - call viscous_coeffs_and_rep_pressure (strength(i,j), DminTarea(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) - + call visccoeff_replpress (strength(i,j) , DminTarea(i,j) , & + Deltane , zetax2 (i,j,1), & + etax2 (i,j,1), rep_prs (i,j,1), & + capping) + call visccoeff_replpress (strength(i,j) , DminTarea(i,j) , & + Deltanw , zetax2 (i,j,2), & + etax2 (i,j,2), rep_prs (i,j,2), & + capping) + call visccoeff_replpress (strength(i,j) , DminTarea(i,j) , & + Deltasw , zetax2 (i,j,3), & + etax2 (i,j,3), rep_prs (i,j,3), & + capping) + call visccoeff_replpress (strength(i,j) , DminTarea(i,j) , & + Deltase , zetax2 (i,j,4), & + etax2 (i,j,4), rep_prs (i,j,4), & + capping) + !----------------------------------------------------------------- ! the stresses ! kg/s^2 ! (1) northeast, (2) northwest, (3) southwest, (4) southeast diff --git a/cicecore/cicedynB/general/ice_init.F90 b/cicecore/cicedynB/general/ice_init.F90 index 4c7817a3d..da6478069 100644 --- a/cicecore/cicedynB/general/ice_init.F90 +++ b/cicecore/cicedynB/general/ice_init.F90 @@ -111,7 +111,8 @@ subroutine input_data k1, k2, alphab, threshold_hw, Ktens, & e_yieldcurve, e_plasticpot, coriolis, & ssh_stress, kridge, brlx, arlx, & - deltaminEVP, deltaminVP, capping + deltaminEVP, deltaminVP, capping, & + elasticDamp use ice_dyn_vp, only: maxits_nonlin, precond, dim_fgmres, dim_pgmres, maxits_fgmres, & maxits_pgmres, monitor_nonlin, monitor_fgmres, & @@ -215,7 +216,7 @@ subroutine input_data namelist /dynamics_nml/ & kdyn, ndte, revised_evp, yield_curve, & - evp_algorithm, & + evp_algorithm, elasticDamp, & brlx, arlx, ssh_stress, & advection, coriolis, kridge, ktransport, & kstrength, krdg_partic, krdg_redist, mu_rdg, & @@ -357,11 +358,12 @@ subroutine input_data ndtd = 1 ! dynamic time steps per thermodynamic time step ndte = 120 ! subcycles per dynamics timestep: ndte=dt_dyn/dte evp_algorithm = 'standard_2d' ! EVP kernel (=standard_2d: standard cice evp; =shared_mem_1d: 1d shared memory and no mpi. if more mpi processors then executed on master + elasticDamp = 0.36_dbl_kind ! coefficient for calculating the parameter E pgl_global_ext = .false. ! if true, init primary grid lengths (global ext.) brlx = 300.0_dbl_kind ! revised_evp values. Otherwise overwritten in ice_dyn_shared arlx = 300.0_dbl_kind ! revised_evp values. Otherwise overwritten in ice_dyn_shared revised_evp = .false. ! if true, use revised procedure for evp dynamics - yield_curve = 'ellipse' + yield_curve = 'ellipse' ! yield curve kstrength = 1 ! 1 = Rothrock 75 strength, 0 = Hibler 79 Pstar = 2.75e4_dbl_kind ! constant in Hibler strength formula (kstrength = 0) Cstar = 20._dbl_kind ! constant in Hibler strength formula (kstrength = 0) @@ -833,6 +835,7 @@ subroutine input_data call broadcast_scalar(ndtd, master_task) call broadcast_scalar(ndte, master_task) call broadcast_scalar(evp_algorithm, master_task) + call broadcast_scalar(elasticDamp, master_task) call broadcast_scalar(pgl_global_ext, master_task) call broadcast_scalar(brlx, master_task) call broadcast_scalar(arlx, master_task) @@ -1666,10 +1669,10 @@ subroutine input_data endif if (kdyn == 1 .or. kdyn == 3) then - write(nu_diag,1030) ' yield_curve = ', trim(yield_curve) + write(nu_diag,1030) ' yield_curve = ', trim(yield_curve), ' : yield curve' if (trim(yield_curve) == 'ellipse') & - write(nu_diag,1002) ' e_yieldcurve = ', e_yieldcurve, ' : aspect ratio of yield curve' - write(nu_diag,1002) ' e_plasticpot = ', e_plasticpot, ' : aspect ratio of plastic potential' + write(nu_diag,1002) ' e_yieldcurve = ', e_yieldcurve, ' : aspect ratio of yield curve' + write(nu_diag,1002) ' e_plasticpot = ', e_plasticpot, ' : aspect ratio of plastic potential' endif if (kdyn == 1) then @@ -1680,6 +1683,8 @@ subroutine input_data write(nu_diag,1002) ' capping = ', capping, ' : capping method for viscous coefficients' endif + write(nu_diag,1002) ' elasticDamp = ', elasticDamp, ' : coefficient for calculating the parameter E' + if (trim(coriolis) == 'latitude') then tmpstr2 = ' : latitude-dependent Coriolis parameter' elseif (trim(coriolis) == 'contant') then @@ -1729,7 +1734,7 @@ subroutine input_data endif endif if (grid_ice == 'C' .or. grid_ice == 'CD') then - write(nu_diag,1030) 'viscous coeff method (U point) = ', trim(visc_coeff_method) + write(nu_diag,1030) ' visc_coeff_method= ', trim(visc_coeff_method),' : viscous coeff method (U point)' endif write(nu_diag,1002) ' Ktens = ', Ktens, ' : tensile strength factor' @@ -2341,7 +2346,7 @@ subroutine input_data 1000 format (a20,1x,f13.6,1x,a) ! float 1002 format (a20,5x,f9.2,1x,a) - 1003 format (a20,1x,G11.4,1x,a) + 1003 format (a20,1x,G13.4,1x,a) 1009 format (a20,1x,d13.6,1x,a) 1010 format (a20,8x,l6,1x,a) ! logical 1011 format (a20,1x,l6) diff --git a/configuration/scripts/tests/gridsys_suite.ts b/configuration/scripts/tests/gridsys_suite.ts index 01132ff83..19f512f14 100644 --- a/configuration/scripts/tests/gridsys_suite.ts +++ b/configuration/scripts/tests/gridsys_suite.ts @@ -9,15 +9,15 @@ smoke gbox80 3x3 boxwall smoke gbox80 2x2 boxwallblock smoke gbox80 1x1 boxslotcyl smoke gbox80 2x4 boxnodyn -smoke gbox80 2x2 boxsymn -smoke gbox80 4x2 boxsyme -smoke gbox80 4x1 boxsymne -smoke gbox80 2x2 boxsymn,kmtislands -smoke gbox80 4x1 boxsyme,kmtislands -smoke gbox80 4x2 boxsymne,kmtislands -smoke gbox80 8x1 boxislandsn -smoke gbox80 4x2 boxislandse -smoke gbox80 2x4 boxislandsne +smoke gbox80 2x2 boxsymn,run1day +smoke gbox80 4x2 boxsyme,run1day +smoke gbox80 4x1 boxsymne,run1day +smoke gbox80 2x2 boxsymn,run1day,kmtislands +smoke gbox80 4x1 boxsyme,run1day,kmtislands +smoke gbox80 4x2 boxsymne,run1day,kmtislands +smoke gbox80 8x1 boxislandsn,run1day +smoke gbox80 4x2 boxislandse,run1day +smoke gbox80 2x4 boxislandsne,run1day smoke gx3 1x1x100x116x1 reprosum,run10day smoke gx3 1x1x25x29x16 reprosum,run10day,dwblockall smoke_gx3_1x1x100x116x1_reprosum_run10day smoke gx3 1x1x5x4x580 reprosum,run10day,dwblockall smoke_gx3_1x1x100x116x1_reprosum_run10day @@ -37,15 +37,15 @@ smoke gbox80 3x3 boxwall,gridcd smoke gbox80 2x2 boxwallblock,gridcd smoke gbox80 1x1 boxslotcyl,gridcd smoke gbox80 2x4 boxnodyn,gridcd -smoke gbox80 2x2 boxsymn,gridcd -smoke gbox80 4x2 boxsyme,gridcd -smoke gbox80 4x1 boxsymne,gridcd -smoke gbox80 2x2 boxsymn,kmtislands,gridcd -smoke gbox80 4x1 boxsyme,kmtislands,gridcd -smoke gbox80 4x2 boxsymne,kmtislands,gridcd -smoke gbox80 8x1 boxislandsn,gridcd -smoke gbox80 4x2 boxislandse,gridcd -smoke gbox80 2x4 boxislandsne,gridcd +smoke gbox80 2x2 boxsymn,run1day,gridcd +smoke gbox80 4x2 boxsyme,run1day,gridcd +smoke gbox80 4x1 boxsymne,run1day,gridcd +smoke gbox80 2x2 boxsymn,run1day,kmtislands,gridcd +smoke gbox80 4x1 boxsyme,run1day,kmtislands,gridcd +smoke gbox80 4x2 boxsymne,run1day,kmtislands,gridcd +smoke gbox80 8x1 boxislandsn,run1day,gridcd +smoke gbox80 4x2 boxislandse,run1day,gridcd +smoke gbox80 2x4 boxislandsne,run1day,gridcd smoke gx3 1x1x100x116x1 reprosum,run10day,gridcd smoke gx3 1x1x25x29x16 reprosum,run10day,dwblockall,gridcd smoke_gx3_1x1x100x116x1_gridcd_reprosum_run10day smoke gx3 1x1x5x4x580 reprosum,run10day,dwblockall,gridcd smoke_gx3_1x1x100x116x1_gridcd_reprosum_run10day @@ -65,15 +65,15 @@ smoke gbox80 3x3 boxwall,gridc smoke gbox80 2x2 boxwallblock,gridc smoke gbox80 1x1 boxslotcyl,gridc smoke gbox80 2x4 boxnodyn,gridc -smoke gbox80 2x2 boxsymn,gridc -smoke gbox80 4x2 boxsyme,gridc -smoke gbox80 4x1 boxsymne,gridc -smoke gbox80 2x2 boxsymn,kmtislands,gridc -smoke gbox80 4x1 boxsyme,kmtislands,gridc -smoke gbox80 4x2 boxsymne,kmtislands,gridc -smoke gbox80 8x1 boxislandsn,gridc -smoke gbox80 4x2 boxislandse,gridc -smoke gbox80 2x4 boxislandsne,gridc +smoke gbox80 2x2 boxsymn,run1day,gridc +smoke gbox80 4x2 boxsyme,run1day,gridc +smoke gbox80 4x1 boxsymne,run1day,gridc +smoke gbox80 2x2 boxsymn,run1day,kmtislands,gridc +smoke gbox80 4x1 boxsyme,run1day,kmtislands,gridc +smoke gbox80 4x2 boxsymne,run1day,kmtislands,gridc +smoke gbox80 8x1 boxislandsn,run1day,gridc +smoke gbox80 4x2 boxislandse,run1day,gridc +smoke gbox80 2x4 boxislandsne,run1day,gridc smoke gx3 1x1x100x116x1 reprosum,run10day,gridc smoke gx3 1x1x25x29x16 reprosum,run10day,dwblockall,gridc smoke_gx3_1x1x100x116x1_gridc_reprosum_run10day smoke gx3 1x1x5x4x580 reprosum,run10day,dwblockall,gridc smoke_gx3_1x1x100x116x1_gridc_reprosum_run10day diff --git a/doc/source/cice_index.rst b/doc/source/cice_index.rst index 8ec9c8f4a..b058d1503 100644 --- a/doc/source/cice_index.rst +++ b/doc/source/cice_index.rst @@ -215,7 +215,7 @@ either Celsius or Kelvin units). "etax2", "2 x eta (shear viscous coefficient)", "kg/s" "evap", "evaporative water flux", "kg/m\ :math:`^2`/s" "ew_boundary_type", "type of east-west boundary condition", "" - "eyc", "coefficient for calculating the parameter E, 0\ :math:`<` eyc :math:`<`\ 1", "0.36" + "elasticDamp", "coefficient for calculating the parameter E, 0\ :math:`<` elasticDamp :math:`<`\ 1", "0.36" "e_yieldcurve", "yield curve minor/major axis ratio", "2" "e_plasticpot", "plastic potential minor/major axis ratio", "2" "**F**", "", "" diff --git a/doc/source/science_guide/sg_dynamics.rst b/doc/source/science_guide/sg_dynamics.rst index 9c529b8ec..b627f896d 100644 --- a/doc/source/science_guide/sg_dynamics.rst +++ b/doc/source/science_guide/sg_dynamics.rst @@ -457,7 +457,7 @@ for elastic waves, :math:`\Delta t_e < T < \Delta t`, as .. math:: E = {\zeta\over T}, -where :math:`T=E_\circ\Delta t` and :math:`E_\circ` (eyc) is a tunable +where :math:`T=E_\circ\Delta t` and :math:`E_\circ` (elasticDamp) is a tunable parameter less than one. Including the modification proposed by :cite:`Bouillon13` for equations :eq:`sig2` and :eq:`sig12` in order to improve numerical convergence, the stress equations become .. math:: diff --git a/doc/source/user_guide/ug_case_settings.rst b/doc/source/user_guide/ug_case_settings.rst index 51437ae1e..fb3806bfd 100644 --- a/doc/source/user_guide/ug_case_settings.rst +++ b/doc/source/user_guide/ug_case_settings.rst @@ -424,6 +424,7 @@ dynamics_nml "``dim_pgmres``", "integer", "maximum number of Arnoldi iterations for PGMRES preconditioner", "5" "``e_plasticpot``", "real", "aspect ratio of elliptical plastic potential", "2.0" "``e_yieldcurve``", "real", "aspect ratio of elliptical yield curve", "2.0" + "``elasticDamp``", "real", "elastic damping parameter", "0.36" "``evp_algorithm``", "``standard_2d``", "standard 2d EVP memory parallel solver", "standard_2d" "", "``shared_mem_1d``", "1d shared memory solver", "" "``kdyn``", "``-1``", "dynamics algorithm OFF", "1" diff --git a/doc/source/user_guide/ug_implementation.rst b/doc/source/user_guide/ug_implementation.rst index 624d135c3..9fc5069d1 100644 --- a/doc/source/user_guide/ug_implementation.rst +++ b/doc/source/user_guide/ug_implementation.rst @@ -1026,9 +1026,9 @@ t_e`) is thus .. math:: dte = dt\_dyn/ndte. -A second parameter, :math:`E_\circ` (``eyc``), defines the elastic wave +A second parameter, :math:`E_\circ` (``elasticDamp``), defines the elastic wave damping timescale :math:`T`, described in Section :ref:`dynam`, as -``eyc * dt_dyn``. The forcing terms are not updated during the subcycling. +``elasticDamp * dt_dyn``. The forcing terms are not updated during the subcycling. Given the small step (``dte``) at which the EVP dynamics model is subcycled, the elastic parameter :math:`E` is also limited by stability constraints, as discussed in :cite:`Hunke97`. Linear stability