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

seabed stress - remove if statements #673

Merged
merged 9 commits into from
Feb 23, 2022
Merged

seabed stress - remove if statements #673

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
5c036cc
refactor seabed_stress. Bit for bit
TillRasmussen Dec 7, 2021
4ab3476
Removed if statement from stepu. Results are binary identical, howeve…
TillRasmussen Dec 12, 2021
0122dd3
changed capping from logical to numeric in order to remove if stateme…
TillRasmussen Dec 20, 2021
34b074c
clean dyn_finish, correct intent(inout) to intent(in) for Cw, resse C…
TillRasmussen Dec 22, 2021
dbbad0d
Merge branch 'main' into neweap
TillRasmussen Feb 19, 2022
22893ca
Reolve conflicts after updating main
TillRasmussen Feb 19, 2022
c833f90
modified environment for Freya to accomodate for additional OMP commands
TillRasmussen Feb 21, 2022
7376807
Requested changes after review. Only changed in seabed stress and not…
TillRasmussen Feb 21, 2022
e312c5b
move deformation to subcycling
TillRasmussen Feb 22, 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
35 changes: 16 additions & 19 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 @@ -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
Open in desktop
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
Open in desktop
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