Skip to content

Commit

Permalink
Refactor viscous_coeffs_and_rep_pressure, rename ecy to elasticDamp a…
Browse files Browse the repository at this point in the history 
…nd make namelist (CICE-Consortium#67)

* - Rename ecy to elasticDamp and add to dynamics namelist
- Update gridsys tests, run sym tests 1 day
- Fix minor formatting issues with namelist output in log file

* update documentation for elasticDamp

* Clean up viscous_coeffs_and_rep_pressure, improve code reuse, rename subroutines
Covert grid_ice "select case" to "if"

* clean up arg list
  • Loading branch information
@apcraig
apcraig authored Mar 15, 2022
1 parent 762b655 commit 7a2d071
Show file tree
Hide file tree
Showing 9 changed files with 246 additions and 330 deletions.
111 changes: 48 additions & 63 deletions cicecore/cicedynB/dynamics/ice_dyn_evp.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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)
Expand Down Expand Up @@ -649,7 +648,7 @@ subroutine evp (dt)
!$OMP END PARALLEL DO
endif

end select
endif

endif

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)'

Expand Down Expand Up @@ -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

!-----------------------------------------------------------------
Expand Down Expand Up @@ -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) ) &
Expand All @@ -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) ) &
Expand All @@ -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

Expand Down
Loading

0 comments on commit 7a2d071

Please sign in to comment.