Skip to content

Commit

Permalink
Merge pull request #21 from apcraig/evp1d_ghteghtn
Browse files Browse the repository at this point in the history 
Update evp1d implementation
  • Loading branch information
@TillRasmussen
TillRasmussen authored Nov 15, 2023
2 parents 191fd2c + 7d50433 commit 7f3703c
Show file tree
Hide file tree
Showing 21 changed files with 377 additions and 370 deletions.
16 changes: 8 additions & 8 deletions cicecore/cicedyn/dynamics/ice_dyn_core1d.F90
View file
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
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

0 comments on commit 7f3703c

Please sign in to comment.