Skip to content

Commit

Permalink
seabed stress - remove if statements (#673)
Browse files Browse the repository at this point in the history 
* refactor seabed_stress. Bit for bit

* Removed if statement from stepu. Results are binary identical, however taubx and tauby is updated on all iterations instead of just the last one. Not used within iteration

* changed capping from logical to numeric in order to remove if statement. Moved call to deformation out of loop

* clean dyn_finish, correct intent(inout) to intent(in) for Cw, resse Cb in stepu, remove if from seabed_stress_LKD

* Reolve conflicts after updating main

* modified environment for Freya to accomodate for additional OMP commands

* Requested changes after review. Only changed in seabed stress and not bit for bit if cor=0.0

added legacy comment in ice_dyn_finish

* move deformation to subcycling
  • Loading branch information
@TillRasmussen
TillRasmussen authored Feb 23, 2022
1 parent bc210be commit a59fa47
Show file tree
Hide file tree
Showing 7 changed files with 142 additions and 174 deletions.
35 changes: 16 additions & 19 deletions cicecore/cicedynB/dynamics/ice_dyn_eap.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -399,29 +399,29 @@ subroutine eap (dt)
!-----------------------------------------------------------------

if (seabed_stress) then
if ( seabed_stress_method == 'LKD' ) then
!$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime)
do iblk = 1, nblocks
call seabed_stress_factor_LKD (nx_block, ny_block, &
icellu (iblk), &
indxui(:,iblk), indxuj(:,iblk), &
vice(:,:,iblk), aice(:,:,iblk), &
hwater(:,:,iblk), Tbu(:,:,iblk))
enddo
!$OMP END PARALLEL DO

!$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime)
do iblk = 1, nblocks

if ( seabed_stress_method == 'LKD' ) then

call seabed_stress_factor_LKD (nx_block, ny_block, &
icellu (iblk), &
indxui(:,iblk), indxuj(:,iblk), &
vice(:,:,iblk), aice(:,:,iblk), &
hwater(:,:,iblk), Tbu(:,:,iblk))

elseif ( seabed_stress_method == 'probabilistic' ) then
elseif ( seabed_stress_method == 'probabilistic' ) then
!$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime)
do iblk = 1, nblocks

call seabed_stress_factor_prob (nx_block, ny_block, &
icellt(iblk), indxti(:,iblk), indxtj(:,iblk), &
icellu(iblk), indxui(:,iblk), indxuj(:,iblk), &
aicen(:,:,:,iblk), vicen(:,:,:,iblk), &
hwater(:,:,iblk), Tbu(:,:,iblk))
endif

enddo
!$OMP END PARALLEL DO
enddo
!$OMP END PARALLEL DO
endif
endif

do ksub = 1,ndte ! subcycling
Expand Down Expand Up @@ -476,7 +476,6 @@ subroutine eap (dt)
call stepu (nx_block, ny_block, &
icellu (iblk), Cdn_ocn (:,:,iblk), &
indxui (:,iblk), indxuj (:,iblk), &
ksub, &
aiu (:,:,iblk), strtmp (:,:,:), &
uocn (:,:,iblk), vocn (:,:,iblk), &
waterx (:,:,iblk), watery (:,:,iblk), &
Expand Down Expand Up @@ -546,8 +545,6 @@ subroutine eap (dt)
uvel (:,:,iblk), vvel (:,:,iblk), &
uocn (:,:,iblk), vocn (:,:,iblk), &
aiu (:,:,iblk), fm (:,:,iblk), &
strintx (:,:,iblk), strinty (:,:,iblk), &
strairx (:,:,iblk), strairy (:,:,iblk), &
strocnx (:,:,iblk), strocny (:,:,iblk), &
strocnxT(:,:,iblk), strocnyT(:,:,iblk))

Expand Down
119 changes: 51 additions & 68 deletions cicecore/cicedynB/dynamics/ice_dyn_evp.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -98,7 +98,7 @@ subroutine evp (dt)
use ice_dyn_evp_1d, only: ice_dyn_evp_1d_copyin, ice_dyn_evp_1d_kernel, &
ice_dyn_evp_1d_copyout
use ice_dyn_shared, only: evp_algorithm, stack_velocity_field, unstack_velocity_field

use ice_dyn_shared, only: deformations
real (kind=dbl_kind), intent(in) :: &
dt ! time step

Expand Down Expand Up @@ -330,29 +330,29 @@ subroutine evp (dt)
!-----------------------------------------------------------------

if (seabed_stress) then

!$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime)
do iblk = 1, nblocks

if ( seabed_stress_method == 'LKD' ) then

if ( seabed_stress_method == 'LKD' ) then
!$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime)
do iblk = 1, nblocks
call seabed_stress_factor_LKD (nx_block, ny_block, &
icellu (iblk), &
indxui(:,iblk), indxuj(:,iblk), &
vice(:,:,iblk), aice(:,:,iblk), &
hwater(:,:,iblk), Tbu(:,:,iblk))
enddo
!$OMP END PARALLEL DO

elseif ( seabed_stress_method == 'probabilistic' ) then

elseif ( seabed_stress_method == 'probabilistic' ) then
!$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime)
do iblk = 1, nblocks
call seabed_stress_factor_prob (nx_block, ny_block, &
icellt(iblk), indxti(:,iblk), indxtj(:,iblk), &
icellu(iblk), indxui(:,iblk), indxuj(:,iblk), &
aicen(:,:,:,iblk), vicen(:,:,:,iblk), &
hwater(:,:,iblk), Tbu(:,:,iblk))
endif

enddo
!$OMP END PARALLEL DO
enddo
!$OMP END PARALLEL DO

endif
endif

call ice_timer_start(timer_evp_2d)
Expand Down Expand Up @@ -401,35 +401,48 @@ subroutine evp (dt)
do iblk = 1, nblocks

! if (trim(yield_curve) == 'ellipse') then
call stress (nx_block, ny_block, &
ksub, icellt(iblk), &
indxti (:,iblk), indxtj (:,iblk), &
uvel (:,:,iblk), vvel (:,:,iblk), &
dxt (:,:,iblk), dyt (:,:,iblk), &
dxhy (:,:,iblk), dyhx (:,:,iblk), &
cxp (:,:,iblk), cyp (:,:,iblk), &
cxm (:,:,iblk), cym (:,:,iblk), &
tarear (:,:,iblk), tinyarea (:,:,iblk), &
strength (:,:,iblk), &
stressp_1 (:,:,iblk), stressp_2 (:,:,iblk), &
stressp_3 (:,:,iblk), stressp_4 (:,:,iblk), &
stressm_1 (:,:,iblk), stressm_2 (:,:,iblk), &
stressm_3 (:,:,iblk), stressm_4 (:,:,iblk), &
stress12_1(:,:,iblk), stress12_2(:,:,iblk), &
stress12_3(:,:,iblk), stress12_4(:,:,iblk), &
shear (:,:,iblk), divu (:,:,iblk), &
rdg_conv (:,:,iblk), rdg_shear (:,:,iblk), &
strtmp (:,:,:) )
! endif ! yield_curve
call stress (nx_block, ny_block, &
icellt(iblk), &
indxti (:,iblk), indxtj (:,iblk), &
uvel (:,:,iblk), vvel (:,:,iblk), &
dxt (:,:,iblk), dyt (:,:,iblk), &
dxhy (:,:,iblk), dyhx (:,:,iblk), &
cxp (:,:,iblk), cyp (:,:,iblk), &
cxm (:,:,iblk), cym (:,:,iblk), &
tinyarea (:,:,iblk), &
strength (:,:,iblk), &
stressp_1 (:,:,iblk), stressp_2 (:,:,iblk), &
stressp_3 (:,:,iblk), stressp_4 (:,:,iblk), &
stressm_1 (:,:,iblk), stressm_2 (:,:,iblk), &
stressm_3 (:,:,iblk), stressm_4 (:,:,iblk), &
stress12_1(:,:,iblk), stress12_2(:,:,iblk), &
stress12_3(:,:,iblk), stress12_4(:,:,iblk), &
strtmp (:,:,:) )
! endif

!-----------------------------------------------------------------
! on last subcycle, save quantities for mechanical redistribution
!-----------------------------------------------------------------
if (ksub == ndte) then
call deformations (nx_block, ny_block , &
icellt(iblk) , &
indxti(:,iblk) , indxtj(:,iblk) , &
uvel(:,:,iblk) , vvel(:,:,iblk) , &
dxt(:,:,iblk) , dyt(:,:,iblk) , &
cxp(:,:,iblk) , cyp(:,:,iblk) , &
cxm(:,:,iblk) , cym(:,:,iblk) , &
tarear(:,:,iblk) , &
shear(:,:,iblk) , divu(:,:,iblk) , &
rdg_conv(:,:,iblk), rdg_shear(:,:,iblk) )
endif

!-----------------------------------------------------------------
! momentum equation
!-----------------------------------------------------------------

call stepu (nx_block, ny_block, &
icellu (iblk), Cdn_ocn (:,:,iblk), &
indxui (:,iblk), indxuj (:,iblk), &
ksub, &
aiu (:,:,iblk), strtmp (:,:,:), &
uocn (:,:,iblk), vocn (:,:,iblk), &
waterx (:,:,iblk), watery (:,:,iblk), &
Expand Down Expand Up @@ -547,8 +560,6 @@ subroutine evp (dt)
uvel (:,:,iblk), vvel (:,:,iblk), &
uocn (:,:,iblk), vocn (:,:,iblk), &
aiu (:,:,iblk), fm (:,:,iblk), &
strintx (:,:,iblk), strinty (:,:,iblk), &
strairx (:,:,iblk), strairy (:,:,iblk), &
strocnx (:,:,iblk), strocny (:,:,iblk), &
strocnxT(:,:,iblk), strocnyT(:,:,iblk))

Expand All @@ -571,30 +582,27 @@ end subroutine evp
! author: Elizabeth C. Hunke, LANL

subroutine stress (nx_block, ny_block, &
ksub, icellt, &
icellt, &
indxti, indxtj, &
uvel, vvel, &
dxt, dyt, &
dxhy, dyhx, &
cxp, cyp, &
cxm, cym, &
tarear, tinyarea, &
tinyarea, &
strength, &
stressp_1, stressp_2, &
stressp_3, stressp_4, &
stressm_1, stressm_2, &
stressm_3, stressm_4, &
stress12_1, stress12_2, &
stress12_3, stress12_4, &
shear, divu, &
rdg_conv, rdg_shear, &
str )

use ice_dyn_shared, only: strain_rates, deformations, viscous_coeffs_and_rep_pressure

integer (kind=int_kind), intent(in) :: &
nx_block, ny_block, & ! block dimensions
ksub , & ! subcycling step
icellt ! no. of cells where icetmask = 1

integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
Expand All @@ -613,20 +621,13 @@ subroutine stress (nx_block, ny_block, &
cxp , & ! 1.5*HTN - 0.5*HTS
cym , & ! 0.5*HTE - 1.5*HTW
cxm , & ! 0.5*HTN - 1.5*HTS
tarear , & ! 1/tarea
tinyarea ! puny*tarea

real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
stressp_1, stressp_2, stressp_3, stressp_4 , & ! sigma11+sigma22
stressm_1, stressm_2, stressm_3, stressm_4 , & ! sigma11-sigma22
stress12_1,stress12_2,stress12_3,stress12_4 ! sigma12

real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
shear , & ! strain rate II component (1/s)
divu , & ! strain rate I component, velocity divergence (1/s)
rdg_conv , & ! convergence term for ridging (1/s)
rdg_shear ! shear term for ridging (1/s)

real (kind=dbl_kind), dimension(nx_block,ny_block,8), intent(out) :: &
str ! stress combinations

Expand Down Expand Up @@ -654,16 +655,15 @@ subroutine stress (nx_block, ny_block, &
str12ew, str12we, str12ns, str12sn , &
strp_tmp, strm_tmp, tmp

logical :: capping ! of the viscous coef
real(kind=dbl_kind),parameter :: capping = c1 ! of the viscous coef

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

!-----------------------------------------------------------------
! Initialize
!-----------------------------------------------------------------

str(:,:,:) = c0
capping = .true. ! could be later included in ice_in

do ij = 1, icellt
i = indxti(ij)
Expand Down Expand Up @@ -869,23 +869,6 @@ subroutine stress (nx_block, ny_block, &

enddo ! ij

!-----------------------------------------------------------------
! on last subcycle, save quantities for mechanical redistribution
!-----------------------------------------------------------------
if (ksub == ndte) then
call deformations (nx_block , ny_block , &
icellt , &
indxti , indxtj , &
uvel , vvel , &
dxt , dyt , &
cxp , cyp , &
cxm , cym , &
tarear , &
shear , divu , &
rdg_conv , rdg_shear )

endif

end subroutine stress

!=======================================================================
Expand Down
6 changes: 2 additions & 4 deletions cicecore/cicedynB/dynamics/ice_dyn_evp_1d.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -869,10 +869,8 @@ subroutine stepu_last(NA_len, rhow, lb, ub, Cw, aiu, uocn, vocn, &
vvel(iw) = (cca * cc2 - ccb * cc1) / ab2

! calculate seabed stress component for outputs
if (seabed_stress) then
taubx(iw) = -uvel(iw) * Tbu(iw) / (sqrt(uold**2 + vold**2) + u0)
tauby(iw) = -vvel(iw) * Tbu(iw) / (sqrt(uold**2 + vold**2) + u0)
end if
taubx(iw) = -uvel(iw) * Tbu(iw) / (sqrt(uold**2 + vold**2) + u0)
tauby(iw) = -vvel(iw) * Tbu(iw) / (sqrt(uold**2 + vold**2) + u0)

end do
#ifdef _OPENACC
Expand Down
Loading

0 comments on commit a59fa47

Please sign in to comment.