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.

Sign up for GitHub

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

allocate c and cd var in evp, reduce number of "if grid_ice". #778

Merged
merged 7 commits into from
Nov 17, 2022
Merged

allocate c and cd var in evp, reduce number of "if grid_ice". #778

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
f12ca31
allocate c and cd var in evp, reduce number of if grid_ice.
TillRasmussen Oct 20, 2022
b16f181
Homogenized call init_dyn in CICE_initMod, rmoved first_call from loo…
TillRasmussen Nov 10, 2022
45b8f01
Merge branch 'main' into evp-refactor
TillRasmussen Nov 10, 2022
f4a41a7
bug fixes, removal of *.loc, moved more fields to be allocatable
TillRasmussen Nov 13, 2022
25f7f03
Merge branch 'evp-refactor' of github.com:TillRasmussen/CICE into evp…
TillRasmussen Nov 13, 2022
780fda9
Allign indents
TillRasmussen Nov 14, 2022
fffe822
read and write when iceumask, icenmask and iceemask are not allocated
TillRasmussen Nov 15, 2022
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
129 changes: 51 additions & 78 deletions cicecore/cicedynB/dynamics/ice_dyn_eap.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@ module ice_dyn_eap
p001, p027, p055, p111, p166, p222, p25, p333
use ice_fileunits, only: nu_diag, nu_dump_eap, nu_restart_eap
use ice_exit, only: abort_ice
use ice_flux, only: rdg_shear
! use ice_timers, only: &
! ice_timer_start, ice_timer_stop, &
! timer_tmp1, timer_tmp2, timer_tmp3, timer_tmp4, &
Expand All @@ -36,8 +37,7 @@ module ice_dyn_eap

implicit none
private
public :: eap, init_eap, write_restart_eap, read_restart_eap, &
alloc_dyn_eap
public :: eap, init_eap, write_restart_eap, read_restart_eap

! Look-up table needed for calculating structure tensor
integer (int_kind), parameter :: &
Expand Down Expand Up @@ -71,42 +71,16 @@ module ice_dyn_eap
real (kind=dbl_kind) :: &
puny, pi, pi2, piq, pih

!=======================================================================

contains

!=======================================================================
! Allocate space for all variables
!
subroutine alloc_dyn_eap
real (kind=dbl_kind), parameter :: &
kfriction = 0.45_dbl_kind

integer (int_kind) :: ierr
real (kind=dbl_kind), save :: &
invdx, invdy, invda, invsin

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

allocate( a11_1 (nx_block,ny_block,max_blocks), &
a11_2 (nx_block,ny_block,max_blocks), &
a11_3 (nx_block,ny_block,max_blocks), &
a11_4 (nx_block,ny_block,max_blocks), &
a12_1 (nx_block,ny_block,max_blocks), &
a12_2 (nx_block,ny_block,max_blocks), &
a12_3 (nx_block,ny_block,max_blocks), &
a12_4 (nx_block,ny_block,max_blocks), &
e11 (nx_block,ny_block,max_blocks), &
e12 (nx_block,ny_block,max_blocks), &
e22 (nx_block,ny_block,max_blocks), &
yieldstress11(nx_block,ny_block,max_blocks), &
yieldstress12(nx_block,ny_block,max_blocks), &
yieldstress22(nx_block,ny_block,max_blocks), &
s11 (nx_block,ny_block,max_blocks), &
s12 (nx_block,ny_block,max_blocks), &
s22 (nx_block,ny_block,max_blocks), &
a11 (nx_block,ny_block,max_blocks), &
a12 (nx_block,ny_block,max_blocks), &
stat=ierr)
if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory')
!=======================================================================

end subroutine alloc_dyn_eap
contains

!=======================================================================
! Elastic-anisotropic-plastic dynamics driver
Expand Down Expand Up @@ -134,7 +108,8 @@ subroutine eap (dt)
dyn_prep1, dyn_prep2, stepu, dyn_finish, &
seabed_stress_factor_LKD, seabed_stress_factor_prob, &
seabed_stress_method, seabed_stress, &
stack_fields, unstack_fields, iceTmask, iceUmask
stack_fields, unstack_fields, iceTmask, iceUmask, &
fld2, fld3, fld4
use ice_flux, only: rdg_conv, strairxT, strairyT, &
strairxU, strairyU, uocn, vocn, ss_tltx, ss_tlty, fmU, &
strtltxU, strtltyU, strocnxU, strocnyU, strintxU, strintyU, taubxU, taubyU, &
Expand Down Expand Up @@ -186,11 +161,6 @@ subroutine eap (dt)
umass , & ! total mass of ice and snow (u grid)
umassdti ! mass of U-cell/dte (kg/m^2 s)

real (kind=dbl_kind), allocatable :: &
fld2(:,:,:,:), & ! temporary for stacking fields for halo update
fld3(:,:,:,:), & ! temporary for stacking fields for halo update
fld4(:,:,:,:) ! temporary for stacking fields for halo update

real (kind=dbl_kind), dimension(nx_block,ny_block,8):: &
strtmp ! stress combinations for momentum equation

Expand All @@ -214,10 +184,6 @@ subroutine eap (dt)
! Initialize
!-----------------------------------------------------------------

allocate(fld2(nx_block,ny_block,2,max_blocks))
allocate(fld3(nx_block,ny_block,3,max_blocks))
allocate(fld4(nx_block,ny_block,4,max_blocks))

! This call is needed only if dt changes during runtime.
! call set_evp_parameters (dt)

Expand All @@ -226,7 +192,7 @@ subroutine eap (dt)
do j = 1, ny_block
do i = 1, nx_block
rdg_conv (i,j,iblk) = c0
! rdg_shear(i,j,iblk) = c0
rdg_shear(i,j,iblk) = c0 ! always zero. Could be moved
divu (i,j,iblk) = c0
shear(i,j,iblk) = c0
e11(i,j,iblk) = c0
Expand Down Expand Up @@ -554,7 +520,6 @@ subroutine eap (dt)

enddo ! subcycling

deallocate(fld2,fld3,fld4)
if (maskhalo_dyn) call ice_HaloDestroy(halo_info_mask)

!-----------------------------------------------------------------
Expand Down Expand Up @@ -588,6 +553,8 @@ subroutine init_eap

use ice_blocks, only: nx_block, ny_block
use ice_domain, only: nblocks
use ice_calendar, only: dt_dyn
use ice_dyn_shared, only: init_dyn_shared

! local variables

Expand All @@ -599,7 +566,7 @@ subroutine init_eap
eps6 = 1.0e-6_dbl_kind

integer (kind=int_kind) :: &
ix, iy, iz, ia
ix, iy, iz, ia, ierr

integer (kind=int_kind), parameter :: &
nz = 100
Expand All @@ -609,6 +576,8 @@ subroutine init_eap
da, dx, dy, dz, &
phi

real (kind=dbl_kind) :: invstressconviso

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

call icepack_query_parameters(puny_out=puny, &
Expand All @@ -619,6 +588,31 @@ subroutine init_eap

phi = pi/c12 ! diamond shaped floe smaller angle (default phi = 30 deg)

call init_dyn_shared(dt_dyn)

allocate( a11_1 (nx_block,ny_block,max_blocks), &
a11_2 (nx_block,ny_block,max_blocks), &
a11_3 (nx_block,ny_block,max_blocks), &
a11_4 (nx_block,ny_block,max_blocks), &
a12_1 (nx_block,ny_block,max_blocks), &
a12_2 (nx_block,ny_block,max_blocks), &
a12_3 (nx_block,ny_block,max_blocks), &
a12_4 (nx_block,ny_block,max_blocks), &
e11 (nx_block,ny_block,max_blocks), &
e12 (nx_block,ny_block,max_blocks), &
e22 (nx_block,ny_block,max_blocks), &
yieldstress11(nx_block,ny_block,max_blocks), &
yieldstress12(nx_block,ny_block,max_blocks), &
yieldstress22(nx_block,ny_block,max_blocks), &
s11 (nx_block,ny_block,max_blocks), &
s12 (nx_block,ny_block,max_blocks), &
s22 (nx_block,ny_block,max_blocks), &
a11 (nx_block,ny_block,max_blocks), &
a12 (nx_block,ny_block,max_blocks), &
stat=ierr)
if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory')


!$OMP PARALLEL DO PRIVATE(iblk,i,j) SCHEDULE(runtime)
do iblk = 1, nblocks
do j = 1, ny_block
Expand All @@ -640,6 +634,7 @@ subroutine init_eap
a12_2 (i,j,iblk) = c0
a12_3 (i,j,iblk) = c0
a12_4 (i,j,iblk) = c0
rdg_shear (i,j,iblk) = c0
enddo ! i
enddo ! j
enddo ! iblk
Expand All @@ -657,6 +652,9 @@ subroutine init_eap
zinit = -pih
dy = pi/real(ny_yield-1,kind=dbl_kind)
yinit = -dy
invdx = c1/dx
invdy = c1/dy
invda = c1/da

do ia=1,na_yield
do ix=1,nx_yield
Expand Down Expand Up @@ -712,6 +710,12 @@ subroutine init_eap
enddo
enddo

! Factor to maintain the same stress as in EVP (see Section 3)
! Can be set to 1 otherwise

invstressconviso = c1/(c1+kfriction*kfriction)
invsin = c1/sin(pi2/c12) * invstressconviso

end subroutine init_eap

!=======================================================================
Expand Down Expand Up @@ -1590,37 +1594,19 @@ subroutine update_stress_rdg (ksub, ndte, divu, tension, &
rotstemp11s, rotstemp12s, rotstemp22s, &
sig11, sig12, sig22, &
sgprm11, sgprm12, sgprm22, &
invstressconviso, &
Angle_denom_gamma, Angle_denom_alpha, &
Tany_1, Tany_2, &
x, y, dx, dy, da, &
dtemp1, dtemp2, atempprime, &
kxw, kyw, kaw

real (kind=dbl_kind), save :: &
invdx, invdy, invda, invsin

logical (kind=log_kind), save :: &
first_call = .true.

real (kind=dbl_kind), parameter :: &
kfriction = 0.45_dbl_kind

! tcraig, temporary, should be moved to namelist
! turns on interpolation in stress_rdg
logical(kind=log_kind), parameter :: &
interpolate_stress_rdg = .false.

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

! Factor to maintain the same stress as in EVP (see Section 3)
! Can be set to 1 otherwise

if (first_call) then
invstressconviso = c1/(c1+kfriction*kfriction)
invsin = c1/sin(pi2/c12) * invstressconviso
endif

! compute eigenvalues, eigenvectors and angles for structure tensor, strain rates

! 1) structure tensor
Expand Down Expand Up @@ -1717,17 +1703,6 @@ subroutine update_stress_rdg (ksub, ndte, divu, tension, &
if (y > pi) y = y - pi
if (y < 0) y = y + pi

! Now calculate updated stress tensor

if (first_call) then
dx = pi/real(nx_yield-1,kind=dbl_kind)
dy = pi/real(ny_yield-1,kind=dbl_kind)
da = p5/real(na_yield-1,kind=dbl_kind)
invdx = c1/dx
invdy = c1/dy
invda = c1/da
endif

if (interpolate_stress_rdg) then

! Interpolated lookup
Expand Down Expand Up @@ -1869,8 +1844,6 @@ subroutine update_stress_rdg (ksub, ndte, divu, tension, &
+ rotstemp22s*dtemp22
endif

first_call = .false.

end subroutine update_stress_rdg

!=======================================================================
Expand Down
Loading