Skip to content

Commit

Permalink
Revp v5 (#242)
Browse files Browse the repository at this point in the history 
* Old sandbox got too much new sand. New branch with new revised evp

* Allowed for more numbers in floats upon writeout
  • Loading branch information
@TillRasmussen @apcraig
TillRasmussen authored and apcraig committed Nov 21, 2018
1 parent a1cc9a5 commit 0681335
Show file tree
Hide file tree
Showing 4 changed files with 51 additions and 63 deletions.
26 changes: 13 additions & 13 deletions cicecore/cicedynB/dynamics/ice_dyn_evp.F90
View file
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 @@ -706,24 +706,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
75 changes: 27 additions & 48 deletions cicecore/cicedynB/dynamics/ice_dyn_shared.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,11 +57,12 @@ module ice_dyn_shared
e_ratio , & ! e = EVP ellipse aspect ratio
ecci , & ! 1/e^2
dtei , & ! 1/dte, where dte is subcycling timestep (1/s)
dte2T , & ! dte/2T
! dte2T , & ! dte/2T
denom1 ! constants for stress equation

real (kind=dbl_kind), public :: & ! Bouillon et al relaxation constants
arlx1i , & ! alpha1 for stressp
arlx , & ! alpha for stressp
arlx1i , & ! (inverse of alpha) for stressp
brlx ! beta for momentum

real (kind=dbl_kind), allocatable, public :: &
Expand Down Expand Up @@ -210,67 +211,45 @@ subroutine set_evp_parameters (dt)

! local variables

real (kind=dbl_kind) :: &
Se , & ! stability parameter for revised EVP
xi , & ! stability parameter for revised EVP
gamma , & ! stability parameter for revised EVP
xmin, ymin , & ! minimum grid length for ocean points, m
dte , & ! subcycling timestep for EVP dynamics, s
ecc , & ! (ratio of major to minor ellipse axes)^2
tdamp2 ! 2*(wave damping time scale T)
!real (kind=dbl_kind) :: &
!dte , & ! subcycling timestep for EVP dynamics, s
!ecc , & ! (ratio of major to minor ellipse axes)^2
!tdamp2 ! 2*(wave damping time scale T)

character(len=*), parameter :: subname = '(set_evp_parameters)'

! elastic time step
dte = dt/real(ndte,kind=dbl_kind) ! s
dtei = c1/dte ! 1/s
!dte = dt/real(ndte,kind=dbl_kind) ! s
!dtei = c1/dte ! 1/s
dtei = real(ndte,kind=dbl_kind)/dt

! major/minor axis length ratio, squared
ecc = e_ratio**2
ecci = c1/ecc ! 1/ecc
!ecc = e_ratio**2
!ecci = c1/ecc ! 1/ecc
ecci = c1/e_ratio**2 ! 1/ecc

! constants for stress equation
tdamp2 = c2*eyc*dt ! s
dte2T = dte/tdamp2 ! ellipse (unitless)

! grid min/max
xmin = global_minval(dxt, distrb_info, tmask)
ymin = global_minval(dyt, distrb_info, tmask)
xmin = min(xmin,ymin) ! min(dxt, dyt)

! revised evp parameters
Se = 0.86_dbl_kind ! Se > 0.5
xi = 5.5e-3_dbl_kind ! Sv/Sc < 1
gamma = p25 * 1.e11_dbl_kind * dt ! rough estimate (P/m~10^5/10^3)
!tdamp2 = c2*eyc*dt ! s
!dte2T = dte/tdamp2 or c1/(c2*eyc*real(ndte,kind=dbl_kind)) ! ellipse (unitless)

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

! classic evp parameters (but modified equations)
! arlx1i = dte2T
! brlx = dt*dtei

denom1 = c1
arlx1i = c1/arlx
else ! Hunke, JCP 2013 with modified stress eq
revp = c0
arlx1i = dte2T
brlx = dt*dtei

! revised evp parameters
! arlx1i = c2*xi/Se ! 1/alpha1
! brlx = c2*Se*xi*gamma/xmin**2 ! beta

!arlx1i = dte2T
!arlx = c1/arlx1i
!brlx = dt*dtei
arlx = c2*eyc*real(ndte,kind=dbl_kind)
arlx1i = c1/arlx
brlx = real(ndte,kind=dbl_kind)
denom1 = c1/(c1+arlx1i)
endif
if (my_task == master_task) then
write (nu_diag,*) 'arlx, brlx', c1/arlx1i, brlx
write (nu_diag,*) 'Se, Sv, xi', &
sqrt(brlx/(arlx1i*gamma))*xmin, &
p5*brlx/gamma*xmin**2, &
p5*xmin*sqrt(brlx*arlx1i/gamma)
endif

denom1 = c1/(c1+arlx1i)
write (nu_diag,*) 'arlx, arlxi, brlx, denom1', &
arlx, arlx1i, brlx, denom1
endif

end subroutine set_evp_parameters

Expand Down
11 changes: 9 additions & 2 deletions cicecore/cicedynB/general/ice_init.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -92,7 +92,7 @@ subroutine input_data
dxrect, dyrect
use ice_dyn_shared, only: ndte, kdyn, revised_evp, yield_curve, &
basalstress, k1, Ktens, e_ratio, coriolis, &
kridge, ktransport
kridge, ktransport, brlx, arlx
use ice_transport_driver, only: advection
use ice_restoring, only: restore_ice
#ifdef CESMCOUPLED
Expand Down Expand Up @@ -166,6 +166,7 @@ subroutine input_data

namelist /dynamics_nml/ &
kdyn, ndte, revised_evp, yield_curve, &
brlx, arlx, &
advection, coriolis, kridge, ktransport, &
kstrength, krdg_partic, krdg_redist, mu_rdg, &
e_ratio, Ktens, Cf, basalstress, &
Expand Down Expand Up @@ -277,6 +278,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
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'
kstrength = 1 ! 1 = Rothrock 75 strength, 0 = Hibler 79
Expand Down Expand Up @@ -536,6 +539,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(arlx, 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 @@ -972,6 +977,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) ' arlx = ', arlx
if (kdyn == 1) &
write(nu_diag,*) ' yield_curve = ', &
trim(yield_curve)
Expand Down Expand Up @@ -1186,7 +1193,7 @@ subroutine input_data
endif

1000 format (a30,2x,f9.2) ! a30 to align formatted, unformatted statements
1005 format (a30,2x,f9.6) ! float
1005 format (a30,2x,f12.6) ! float
1010 format (a30,2x,l6) ! logical
1020 format (a30,2x,i6) ! integer
1021 format (a30,2x,a8,i6) ! char, int
Expand Down
2 changes: 2 additions & 0 deletions configuration/scripts/ice_in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -102,6 +102,8 @@
kdyn = 1
ndte = 240
revised_evp = .false.
brlx = 300.0
arlx = 300.0
advection = 'remap'
kstrength = 1
krdg_partic = 1
Expand Down

0 comments on commit 0681335

Please sign in to comment.