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

Jump to bottom

Update evp1d implementation #21

Merged
merged 2 commits into from
Nov 15, 2023
Merged

Update evp1d implementation #21

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
a6321b4
Update of evp1d implementation
apcraig Nov 14, 2023
7d50433
Clean up trailing blanks
apcraig Nov 14, 2023
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
16 changes: 8 additions & 8 deletions cicecore/cicedyn/dynamics/ice_dyn_core1d.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,7 @@ subroutine stress_1d (ee, ne, se, lb, ub, &

use ice_dyn_shared, only: arlx1i, denom1, revp, &
deltaminEVP, visc_replpress
!
!
implicit none
! arguments ------------------------------------------------------------------
integer (kind=int_kind), intent(in) :: lb,ub
Expand Down Expand Up @@ -274,7 +274,7 @@ subroutine stress_1d (ee, ne, se, lb, ub, &
ssigpw = stressp_2(iw) + stressp_3(iw)
ssigp1 =(stressp_1(iw) + stressp_3(iw))*p055
ssigp2 =(stressp_2(iw) + stressp_4(iw))*p055

ssigmn = stressm_1(iw) + stressm_2(iw)
ssigms = stressm_3(iw) + stressm_4(iw)
ssigme = stressm_1(iw) + stressm_4(iw)
Expand Down Expand Up @@ -396,7 +396,7 @@ subroutine strain_rates_1d (tmp_uvel_cc, tmp_vvel_cc, &

real (kind=dbl_kind), intent(in) :: &
tmp_uvel_ee, tmp_vvel_ee, tmp_uvel_se, tmp_vvel_se, &
tmp_uvel_cc, tmp_vvel_cc, tmp_uvel_ne, tmp_vvel_ne
tmp_uvel_cc, tmp_vvel_cc, tmp_uvel_ne, tmp_vvel_ne

real (kind=dbl_kind), intent(in) :: &
dxT , & ! width of T-cell through the middle (m)
Expand All @@ -405,7 +405,7 @@ subroutine strain_rates_1d (tmp_uvel_cc, tmp_vvel_cc, &
cxp , & ! 1.5*HTN - 0.5*HTS
cym , & ! 0.5*HTE - 1.5*HTW
cxm ! 0.5*HTN - 1.5*HTS

real (kind=dbl_kind), intent(out):: & ! at each corner :
divune, divunw, divuse, divusw , & ! divergence
tensionne, tensionnw, tensionse, tensionsw, & ! tension
Expand Down Expand Up @@ -479,7 +479,7 @@ subroutine stepu_1d (lb , ub , &
uarear , &
uvel_init, vvel_init, &
uvel , vvel , &
str1 , str2 , &
str1 , str2 , &
str3 , str4 , &
str5 , str6 , &
str7 , str8 , &
Expand Down Expand Up @@ -569,7 +569,7 @@ subroutine stepu_1d (lb , ub , &
! revp = 0 for classic evp, 1 for revised evp
cca = (brlx + revp)*umassdti(iw) + vrel * cosw + Cb(iw) ! kg/m^2 s
ccb = fm(iw) + sign(c1,fm(iw)) * vrel * sinw ! kg/m^2 s

ab2 = cca**2 + ccb**2

tmp_str2_nw = str2(nw(iw))
Expand All @@ -578,7 +578,7 @@ subroutine stepu_1d (lb , ub , &
tmp_str6_sse = str6(sse(iw))
tmp_str7_nw = str7(nw(iw))
tmp_str8_sw = str8(sw(iw))

! divergence of the internal stress tensor
strintx = uarear(iw)*(str1(iw)+tmp_str2_nw+tmp_str3_sse+tmp_str4_sw)
strinty = uarear(iw)*(str5(iw)+tmp_str6_sse+tmp_str7_nw+tmp_str8_sw)
Expand All @@ -590,7 +590,7 @@ subroutine stepu_1d (lb , ub , &
+ umassdti(iw)*(brlx*vold + revp*vvel_init(iw))
uvel(iw) = (cca*cc1 + ccb*cc2) / ab2 ! m/s
vvel(iw) = (cca*cc2 - ccb*cc1) / ab2

! calculate seabed stress component for outputs
! only needed on last iteration.
enddo
Expand Down
103 changes: 51 additions & 52 deletions cicecore/cicedyn/dynamics/ice_dyn_evp.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -134,9 +134,7 @@ subroutine init_evp
call init_dyn_shared(dt_dyn)

if (evp_algorithm == "shared_mem_1d" ) then
call dyn_evp1d_init(nx_global , ny_global, nx_block, ny_block, &
max_blocks, nghost , dyT , dxT , &
uarear , tmask , G_HTE , G_HTN)
call dyn_evp1d_init
endif

allocate( uocnU (nx_block,ny_block,max_blocks), & ! i ocean current (m/s)
Expand Down Expand Up @@ -203,6 +201,7 @@ subroutine init_evp

end subroutine init_evp

!=======================================================================
#ifdef CICE_IN_NEMO
! Wind stress is set during this routine from the values supplied
! via NEMO (unless calc_strair is true). These values are supplied
Expand Down Expand Up @@ -254,7 +253,7 @@ subroutine evp (dt)
strain_rates_U, &
iceTmask, iceUmask, iceEmask, iceNmask, &
dyn_haloUpdate, fld2, fld3, fld4
use ice_dyn_evp1d, only: dyn_evp1d_run
use ice_dyn_evp1d, only: dyn_evp1d_run

real (kind=dbl_kind), intent(in) :: &
dt ! time step
Expand Down Expand Up @@ -801,20 +800,20 @@ subroutine evp (dt)

if (grid_ice == "B") then

if (evp_algorithm == "shared_mem_1d" ) then
if (evp_algorithm == "shared_mem_1d" ) then

call dyn_evp1d_run(stressp_1 , stressp_2, stressp_3 , stressp_4 , &
stressm_1 , stressm_2 , stressm_3 , stressm_4 , &
stress12_1, stress12_2, stress12_3, stress12_4, &
strength , &
cdn_ocnU , aiu , uocnU , vocnU , &
waterxU , wateryU , forcexU , forceyU , &
umassdti , fmU , strintxU , strintyU , &
Tbu , taubxU , taubyU , uvel , &
vvel , icetmask , iceUmask)
call dyn_evp1d_run(stressp_1 , stressp_2, stressp_3 , stressp_4 , &
stressm_1 , stressm_2 , stressm_3 , stressm_4 , &
stress12_1, stress12_2, stress12_3, stress12_4, &
strength , &
cdn_ocnU , aiu , uocnU , vocnU , &
waterxU , wateryU , forcexU , forceyU , &
umassdti , fmU , strintxU , strintyU , &
Tbu , taubxU , taubyU , uvel , &
vvel , icetmask , iceUmask)

else ! evp_algorithm == standard_2d (Standard CICE)
do ksub = 1,ndte ! subcycling
else ! evp_algorithm == standard_2d (Standard CICE)
do ksub = 1,ndte ! subcycling

!$OMP PARALLEL DO PRIVATE(iblk,strtmp) SCHEDULE(runtime)
do iblk = 1, nblocks
Expand Down Expand Up @@ -868,26 +867,26 @@ subroutine evp (dt)
uvel, vvel)

enddo ! sub cycling
endif ! evp algorithm
endif ! evp algorithm

!-----------------------------------------------------------------
! save quantities for mechanical redistribution
!-----------------------------------------------------------------
!-----------------------------------------------------------------
! save quantities for mechanical redistribution
!-----------------------------------------------------------------

!$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime)
do iblk = 1, nblocks
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) )
enddo
!$OMP END PARALLEL DO
!$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime)
do iblk = 1, nblocks
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) )
enddo
!$OMP END PARALLEL DO

elseif (grid_ice == "C") then

Expand Down Expand Up @@ -1079,8 +1078,8 @@ subroutine evp (dt)

do ksub = 1,ndte ! subcycling

!$OMP PARALLEL DO PRIVATE(iblk)
do iblk = 1, nblocks
!$OMP PARALLEL DO PRIVATE(iblk)
do iblk = 1, nblocks
call stressCD_T (nx_block , ny_block , &
icellT (iblk), &
indxTi (:,iblk), indxTj (:,iblk), &
Expand All @@ -1094,27 +1093,27 @@ subroutine evp (dt)
stresspT (:,:,iblk), stressmT (:,:,iblk), &
stress12T (:,:,iblk) )

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

! calls ice_haloUpdate, controls bundles and masks
call dyn_haloUpdate (halo_info, halo_info_mask, &
field_loc_center, field_type_scalar, &
zetax2T, etax2T)
! calls ice_haloUpdate, controls bundles and masks
call dyn_haloUpdate (halo_info, halo_info_mask, &
field_loc_center, field_type_scalar, &
zetax2T, etax2T)

if (visc_method == 'avg_strength') then
if (visc_method == 'avg_strength') then
call grid_average_X2Y('S', strength, 'T', strengthU, 'U')
elseif (visc_method == 'avg_zeta') then
elseif (visc_method == 'avg_zeta') then
call grid_average_X2Y('S', zetax2T , 'T', zetax2U , 'U')
call grid_average_X2Y('S', etax2T , 'T', etax2U , 'U')
endif

!$OMP PARALLEL DO PRIVATE(iblk)
do iblk = 1, nblocks
!-----------------------------------------------------------------
! strain rates at U point
! NOTE these are actually strain rates * area (m^2/s)
!-----------------------------------------------------------------
endif

!$OMP PARALLEL DO PRIVATE(iblk)
do iblk = 1, nblocks
!-----------------------------------------------------------------
! strain rates at U point
! NOTE these are actually strain rates * area (m^2/s)
!-----------------------------------------------------------------
call strain_rates_U (nx_block , ny_block , &
icellU (iblk), &
indxUi (:,iblk), indxUj (:,iblk), &
Expand Down
Loading
Loading