Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Revised EVP updated #184

Closed
wants to merge 3 commits into from
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
26 changes: 13 additions & 13 deletions cicecore/cicedynB/dynamics/ice_dyn_evp.F90
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@ module ice_dyn_evp
p2, p222, p25, p333, p5, c1
use ice_dyn_shared, only: stepu, dyn_prep1, dyn_prep2, dyn_finish, &
ndte, yield_curve, ecci, denom1, arlx1i, fcor_blk, uvel_init, &
vvel_init, basal_stress_coeff, basalstress, Ktens
vvel_init, basal_stress_coeff, basalstress, Ktens, revp
use ice_fileunits, only: nu_diag
use ice_exit, only: abort_ice
use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted
Expand Down Expand Up @@ -710,24 +710,24 @@ subroutine stress (nx_block, ny_block, &
! (1) northeast, (2) northwest, (3) southwest, (4) southeast
!-----------------------------------------------------------------

stressp_1(i,j) = (stressp_1(i,j) + c1ne*(divune*(c1+Ktens) - Deltane*(c1-Ktens))) &
stressp_1(i,j) = (stressp_1(i,j)*(c1-arlx1i*revp) + c1ne*(divune*(c1+Ktens) - Deltane*(c1-Ktens))) &
* denom1
stressp_2(i,j) = (stressp_2(i,j) + c1nw*(divunw*(c1+Ktens) - Deltanw*(c1-Ktens))) &
stressp_2(i,j) = (stressp_2(i,j)*(c1-arlx1i*revp) + c1nw*(divunw*(c1+Ktens) - Deltanw*(c1-Ktens))) &
* denom1
stressp_3(i,j) = (stressp_3(i,j) + c1sw*(divusw*(c1+Ktens) - Deltasw*(c1-Ktens))) &
stressp_3(i,j) = (stressp_3(i,j)*(c1-arlx1i*revp) + c1sw*(divusw*(c1+Ktens) - Deltasw*(c1-Ktens))) &
* denom1
stressp_4(i,j) = (stressp_4(i,j) + c1se*(divuse*(c1+Ktens) - Deltase*(c1-Ktens))) &
stressp_4(i,j) = (stressp_4(i,j)*(c1-arlx1i*revp) + c1se*(divuse*(c1+Ktens) - Deltase*(c1-Ktens))) &
* denom1

stressm_1(i,j) = (stressm_1(i,j) + c0ne*tensionne*(c1+Ktens)) * denom1
stressm_2(i,j) = (stressm_2(i,j) + c0nw*tensionnw*(c1+Ktens)) * denom1
stressm_3(i,j) = (stressm_3(i,j) + c0sw*tensionsw*(c1+Ktens)) * denom1
stressm_4(i,j) = (stressm_4(i,j) + c0se*tensionse*(c1+Ktens)) * denom1
stressm_1(i,j) = (stressm_1(i,j)*(c1-arlx1i*revp) + c0ne*tensionne*(c1+Ktens)) * denom1
stressm_2(i,j) = (stressm_2(i,j)*(c1-arlx1i*revp) + c0nw*tensionnw*(c1+Ktens)) * denom1
stressm_3(i,j) = (stressm_3(i,j)*(c1-arlx1i*revp) + c0sw*tensionsw*(c1+Ktens)) * denom1
stressm_4(i,j) = (stressm_4(i,j)*(c1-arlx1i*revp) + c0se*tensionse*(c1+Ktens)) * denom1

stress12_1(i,j) = (stress12_1(i,j) + c0ne*shearne*p5*(c1+Ktens)) * denom1
stress12_2(i,j) = (stress12_2(i,j) + c0nw*shearnw*p5*(c1+Ktens)) * denom1
stress12_3(i,j) = (stress12_3(i,j) + c0sw*shearsw*p5*(c1+Ktens)) * denom1
stress12_4(i,j) = (stress12_4(i,j) + c0se*shearse*p5*(c1+Ktens)) * denom1
stress12_1(i,j) = (stress12_1(i,j)*(c1-arlx1i*revp) + c0ne*shearne*p5*(c1+Ktens)) * denom1
stress12_2(i,j) = (stress12_2(i,j)*(c1-arlx1i*revp) + c0nw*shearnw*p5*(c1+Ktens)) * denom1
stress12_3(i,j) = (stress12_3(i,j)*(c1-arlx1i*revp) + c0sw*shearsw*p5*(c1+Ktens)) * denom1
stress12_4(i,j) = (stress12_4(i,j)*(c1-arlx1i*revp) + c0se*shearse*p5*(c1+Ktens)) * denom1

!-----------------------------------------------------------------
! Eliminate underflows.
Expand Down
11 changes: 7 additions & 4 deletions cicecore/cicedynB/dynamics/ice_dyn_shared.F90
Original file line number Diff line number Diff line change
Expand Up @@ -242,8 +242,9 @@ subroutine set_evp_parameters (dt)

if (revised_evp) then ! Bouillon et al, Ocean Mod 2013
revp = c1
arlx1i = c2*xi/Se ! 1/alpha1
brlx = c2*Se*xi*gamma/xmin**2 ! beta
denom1 = c1
!TAR: According to mail of Martin Loesch: for a ~25km grid:
! brlx=300 , arlx1i = c1/brlx

! classic evp parameters (but modified equations)
! arlx1i = dte2T
Expand All @@ -253,6 +254,8 @@ subroutine set_evp_parameters (dt)
revp = c0
arlx1i = dte2T
brlx = dt*dtei
denom1 = c1/(c1+arlx1i)


! revised evp parameters
! arlx1i = c2*xi/Se ! 1/alpha1
Expand All @@ -267,8 +270,6 @@ subroutine set_evp_parameters (dt)
p5*xmin*sqrt(brlx*arlx1i/gamma)
endif

denom1 = c1/(c1+arlx1i)

end subroutine set_evp_parameters

!=======================================================================
Expand Down Expand Up @@ -526,6 +527,8 @@ subroutine dyn_prep2 (nx_block, ny_block, &
taubx (i,j) = c0
tauby (i,j) = c0

! Don't understand why this differs between revp and none revp.
! However equal nevertheless
if (revp==1) then ! revised evp
stressp_1 (i,j) = c0
stressp_2 (i,j) = c0
Expand Down
9 changes: 8 additions & 1 deletion cicecore/cicedynB/general/ice_init.F90
Original file line number Diff line number Diff line change
Expand Up @@ -89,7 +89,7 @@ subroutine input_data
use ice_grid, only: grid_file, gridcpl_file, kmt_file, grid_type, grid_format, &
dxrect, dyrect
use ice_dyn_shared, only: ndte, kdyn, revised_evp, yield_curve, &
basalstress, Ktens, e_ratio, coriolis, &
basalstress, Ktens, e_ratio, coriolis, brlx, arlx1i, &
kridge, ktransport
use ice_transport_driver, only: advection
use ice_restoring, only: restore_ice
Expand Down Expand Up @@ -164,6 +164,7 @@ subroutine input_data

namelist /dynamics_nml/ &
kdyn, ndte, revised_evp, yield_curve, &
brlx, arlx1i, &
advection, coriolis, kridge, ktransport, &
kstrength, krdg_partic, krdg_redist, mu_rdg, &
e_ratio, Ktens, Cf, basalstress
Expand Down Expand Up @@ -275,6 +276,8 @@ subroutine input_data
kdyn = 1 ! type of dynamics (-1, 0 = off, 1 = evp, 2 = eap)
ndtd = 1 ! dynamic time steps per thermodynamic time step
ndte = 120 ! subcycles per dynamics timestep: ndte=dt_dyn/dte
brlx = 300.0_dbl_kind ! revised_evp values. Otherwise overwritten in ice_dyn_shared
arlx1i = c1/brlx ! revised_evp values. Otherwise overwritten in ice_dyn_shared
revised_evp = .false. ! if true, use revised procedure for evp dynamics
yield_curve = 'ellipse'
kstrength = 1 ! 1 = Rothrock 75 strength, 0 = Hibler 79
Expand Down Expand Up @@ -542,6 +545,8 @@ subroutine input_data
call broadcast_scalar(kdyn, master_task)
call broadcast_scalar(ndtd, master_task)
call broadcast_scalar(ndte, master_task)
call broadcast_scalar(brlx, master_task)
call broadcast_scalar(arlx1i, master_task)
call broadcast_scalar(revised_evp, master_task)
call broadcast_scalar(yield_curve, master_task)
call broadcast_scalar(kstrength, master_task)
Expand Down Expand Up @@ -983,6 +988,8 @@ subroutine input_data
write(nu_diag,1020) ' ndte = ', ndte
write(nu_diag,1010) ' revised_evp = ', &
revised_evp
write(nu_diag,1005) ' brlx = ', brlx
write(nu_diag,1005) ' arlx1i = ', arlx1i
if (kdyn == 1) &
write(nu_diag,*) ' yield_curve = ', &
trim(yield_curve)
Expand Down
2 changes: 2 additions & 0 deletions configuration/scripts/ice_in
Original file line number Diff line number Diff line change
Expand Up @@ -108,6 +108,8 @@
&dynamics_nml
kdyn = 1
ndte = 120
brlx = 300.0
arlx1i = 1.0/300.0
revised_evp = .false.
advection = 'remap'
kstrength = 1
Expand Down