Skip to content

Commit

Permalink
allocate c and cd var in evp, reduce number of "if grid_ice". (#778)
Browse files Browse the repository at this point in the history 
* allocate c and cd var in evp, reduce number of if grid_ice.

* Homogenized call init_dyn in CICE_initMod, rmoved first_call from loop in eap

* bug fixes, removal of *.loc, moved more fields to be allocatable

* Allign indents

* read and write when iceumask, icenmask and iceemask are not allocated
  • Loading branch information
@TillRasmussen
TillRasmussen authored Nov 17, 2022
1 parent 9808b51 commit 99daf69
Show file tree
Hide file tree
Showing 15 changed files with 573 additions and 503 deletions.
129 changes: 51 additions & 78 deletions cicecore/cicedynB/dynamics/ice_dyn_eap.F90
View file
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

0 comments on commit 99daf69

Please sign in to comment.