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.

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
Prev Previous commit
Next Next commit
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
  • Loading branch information
TillRasmussen committed Feb 22, 2022
commit 7376807288e776f64e25aa8a87058b978203be2e
20 changes: 10 additions & 10 deletions cicecore/cicedynB/dynamics/ice_dyn_evp.F90
Original file line number Diff line number Diff line change
Expand Up @@ -483,16 +483,16 @@ subroutine evp (dt)
! on last subcycle, save quantities for mechanical redistribution
!-----------------------------------------------------------------

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) )
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) )


!-----------------------------------------------------------------
Expand Down
18 changes: 14 additions & 4 deletions cicecore/cicedynB/dynamics/ice_dyn_shared.F90
Original file line number Diff line number Diff line change
Expand Up @@ -176,7 +176,7 @@ subroutine init_dyn (dt)
if (trim(coriolis) == 'constant') then
fcor_blk(i,j,iblk) = 1.46e-4_dbl_kind ! Hibler 1979, N. Hem; 1/s
else if (trim(coriolis) == 'zero') then
fcor_blk(i,j,iblk) = 0.0
fcor_blk(i,j,iblk) = c0
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is the change in this line the entire reason the code is not BFB for cor=0?

Copy link
Contributor Author

@TillRasmussen TillRasmussen Feb 21, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yes that is why I put in the comment. I dont think that I reran a test with coriolis=0.0 after the last commit but just a warning if the standard test complains about bit for bit differences.
All other test have been bit for bit with main after Tony's OMP updates.
.

else
fcor_blk(i,j,iblk) = c2*omega*sin(ULAT(i,j,iblk)) ! 1/s
endif
Expand Down Expand Up @@ -760,6 +760,8 @@ subroutine dyn_finish (nx_block, ny_block, &
uvel, vvel, &
uocn, vocn, &
aiu, fm, &
! strintx, strinty, &
! strairx, strairy, &
strocnx, strocny, &
strocnxT, strocnyT)

Expand All @@ -778,6 +780,10 @@ subroutine dyn_finish (nx_block, ny_block, &
vocn , & ! ocean current, y-direction (m/s)
aiu , & ! ice fraction on u-grid
fm ! Coriolis param. * mass in U-cell (kg/s)
! strintx , & ! divergence of internal ice stress, x (N/m^2)
! strinty , & ! divergence of internal ice stress, y (N/m^2)
! strairx , & ! stress on ice by air, x-direction
! strairy ! stress on ice by air, y-direction

real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
strocnx , & ! ice-ocean stress, x-direction
Expand Down Expand Up @@ -905,8 +911,8 @@ subroutine seabed_stress_factor_LKD (nx_block, ny_block, &

hwu = min(hwater(i,j),hwater(i+1,j),hwater(i,j+1),hwater(i+1,j+1))

! if (hwu < threshold_hw) then
docalc_tbu = max(sign(c1,threshold_hw-hwu),c0)
docalc_tbu = merge(c1,c0,hwu < threshold_hw)


au = max(aice(i,j),aice(i+1,j),aice(i,j+1),aice(i+1,j+1))
hu = max(vice(i,j),vice(i+1,j),vice(i,j+1),vice(i+1,j+1))
Expand Down Expand Up @@ -1395,7 +1401,8 @@ subroutine viscous_coeffs_and_rep_pressure (strength, tinyarea, &
tmpcalcne, tmpcalcnw, tmpcalcsw, tmpcalcse

! NOTE: for comp. efficiency 2 x zeta and 2 x eta are used in the code


! if (trim(yield_curve) == 'ellipse') then
tmpcalcne = capping *(strength/max(Deltane, tinyarea))+ &
(c1-capping)* strength/ (Deltane+ tinyarea)
tmpcalcnw = capping *(strength/max(Deltanw, tinyarea))+ &
Expand All @@ -1420,6 +1427,9 @@ subroutine viscous_coeffs_and_rep_pressure (strength, tinyarea, &
zetax2se = (c1+Ktens)*tmpcalcse ! southeast
rep_prsse = (c1-Ktens)*tmpcalcse*Deltase
etax2se = epp2i*zetax2se
! else

! endif

end subroutine viscous_coeffs_and_rep_pressure

Expand Down