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

New 1d evp solver #895

Merged
merged 13 commits into from
Nov 16, 2023
Merged

New 1d evp solver #895

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
6f6e2d2
New 1d evp solver
TillRasmussen Oct 24, 2023
8f1e52d
Small changes incl timer names and inclued private/publice in ice_dyn…
TillRasmussen Nov 6, 2023
4e42526
fixed bug on gnu debug
TillRasmussen Nov 7, 2023
5e13f64
moved halo update to evp1d, added deallocation, fixed bug
TillRasmussen Nov 10, 2023
6a65fb6
fixed deallocation dyn_evp1d
TillRasmussen Nov 10, 2023
435781b
bugfix deallocate
TillRasmussen Nov 10, 2023
6fc1673
Remove gather strintx and strinty
TillRasmussen Nov 11, 2023
191fd2c
removed 4 test with evp1d and c/cd grid
TillRasmussen Nov 12, 2023
a6321b4
Update of evp1d implementation
apcraig Nov 14, 2023
7d50433
Clean up trailing blanks
apcraig Nov 14, 2023
7f3703c
Merge pull request #21 from apcraig/evp1d_ghteghtn
TillRasmussen Nov 15, 2023
db24921
resolved name conflicts
TillRasmussen Nov 15, 2023
7c22713
1d grid var name change
TillRasmussen Nov 15, 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
671 changes: 671 additions & 0 deletions cicecore/cicedyn/dynamics/ice_dyn_core1d.F90
View file
Open in desktop

Large diffs are not rendered by default.

794 changes: 390 additions & 404 deletions cicecore/cicedyn/dynamics/ice_dyn_evp.F90
View file
Open in desktop

Large diffs are not rendered by default.

1,467 changes: 1,467 additions & 0 deletions cicecore/cicedyn/dynamics/ice_dyn_evp1d.F90
View file
Open in desktop

Large diffs are not rendered by default.

1,921 changes: 0 additions & 1,921 deletions cicecore/cicedyn/dynamics/ice_dyn_evp_1d.F90
View file
Open in desktop

This file was deleted.

15 changes: 7 additions & 8 deletions cicecore/cicedyn/dynamics/ice_dyn_shared.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1451,10 +1451,10 @@ subroutine seabed_stress_factor_prob (nx_block, ny_block, &

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

call icepack_query_parameters(rhow_out=rhow, rhoi_out=rhoi)
call icepack_query_parameters(gravit_out=gravit)
call icepack_query_parameters(pi_out=pi)
call icepack_query_parameters(puny_out=puny)
call icepack_query_parameters(rhow_out=rhow, rhoi_out=rhoi,gravit_out=gravit,pi_out=pi,puny_out=puny)
call icepack_warnings_flush(nu_diag)
if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
file=__FILE__, line=__LINE__)

Tbt=c0

Expand Down Expand Up @@ -2302,16 +2302,15 @@ end subroutine strain_rates_U
! by combining tensile strength and a parameterization for grounded ridges.
! J. Geophys. Res. Oceans, 121, 7354-7368.

subroutine visc_replpress(strength, DminArea, Delta, &
zetax2, etax2, rep_prs, capping)
subroutine visc_replpress(strength, DminArea, Delta, &
zetax2, etax2, rep_prs)

real (kind=dbl_kind), intent(in):: &
strength, & !
DminArea !

real (kind=dbl_kind), intent(in):: &
Delta , & !
capping !
Delta

real (kind=dbl_kind), intent(out):: &
zetax2 , & ! bulk viscosity
Expand Down
12 changes: 4 additions & 8 deletions cicecore/cicedyn/dynamics/ice_dyn_vp.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1221,20 +1221,16 @@ subroutine calc_zeta_dPr (nx_block, ny_block, &

call visc_replpress (strength(i,j) , DminTarea(i,j) , &
Deltane , zetax2 (i,j,1), &
etax2 (i,j,1), rep_prs (i,j,1), &
capping)
etax2 (i,j,1), rep_prs (i,j,1))
call visc_replpress (strength(i,j) , DminTarea(i,j) , &
Deltanw , zetax2 (i,j,2), &
etax2 (i,j,2), rep_prs (i,j,2), &
capping)
etax2 (i,j,2), rep_prs (i,j,2))
call visc_replpress (strength(i,j) , DminTarea(i,j) , &
Deltasw , zetax2 (i,j,3), &
etax2 (i,j,3), rep_prs (i,j,3), &
capping)
etax2 (i,j,3), rep_prs (i,j,3))
call visc_replpress (strength(i,j) , DminTarea(i,j) , &
Deltase , zetax2 (i,j,4), &
etax2 (i,j,4), rep_prs (i,j,4), &
capping)
etax2 (i,j,4), rep_prs (i,j,4))

!-----------------------------------------------------------------
! the stresses ! kg/s^2
Expand Down
49 changes: 25 additions & 24 deletions cicecore/cicedyn/dynamics/ice_transport_remap.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -308,12 +308,12 @@ subroutine init_remap
! regions are adjusted to obtain the desired area.
! If false, edgearea is computed in locate_triangles and passed out.
!
! l_fixed_area = .false. has been the default approach in CICE. It is
! used like this for the B-grid. However, idealized tests with the
! C-grid have shown that l_fixed_area = .false. leads to a checkerboard
! pattern in prognostic fields (e.g. aice). Using l_fixed_area = .true.
! l_fixed_area = .false. has been the default approach in CICE. It is
! used like this for the B-grid. However, idealized tests with the
! C-grid have shown that l_fixed_area = .false. leads to a checkerboard
! pattern in prognostic fields (e.g. aice). Using l_fixed_area = .true.
! eliminates the checkerboard pattern in C-grid simulations.
!
!
!-------------------------------------------------------------------

if (grid_ice == 'CD' .or. grid_ice == 'C') then
Expand Down Expand Up @@ -647,11 +647,12 @@ subroutine horizontal_remap (dt, ntrace, &
endif ! nghost

! tcraig, this OMP loop sometimes fails with cce/14.0.3, compiler bug??
!$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block,n, &
!$OMP edgearea_e,edgearea_n,edge,iflux,jflux, &
!$OMP xp,yp,indxing,indxjng,mflxe,mflxn, &
!$OMP mtflxe,mtflxn,triarea,istop,jstop,l_stop) &
!$OMP SCHEDULE(runtime)
! TILL I can trigger the same with ifort (IFORT) 18.0.0 20170811
!TILL !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block,n, &
!TILL !$OMP edgearea_e,edgearea_n,edge,iflux,jflux, &
!TILL !$OMP xp,yp,indxing,indxjng,mflxe,mflxn, &
!TILL !$OMP mtflxe,mtflxn,triarea,istop,jstop,l_stop) &
!TILL !$OMP SCHEDULE(runtime)
do iblk = 1, nblocks

l_stop = .false.
Expand Down Expand Up @@ -865,7 +866,7 @@ subroutine horizontal_remap (dt, ntrace, &
enddo ! n

enddo ! iblk
!$OMP END PARALLEL DO
!TILL !$OMP END PARALLEL DO

end subroutine horizontal_remap

Expand Down Expand Up @@ -1725,7 +1726,7 @@ subroutine locate_triangles (nx_block, ny_block, &
dpy , & ! y coordinates of departure points at cell corners
dxu , & ! E-W dimension of U-cell (m)
dyu , & ! N-S dimension of U-cell (m)
earea , & ! area of E-cell
earea , & ! area of E-cell
narea ! area of N-cell

real (kind=dbl_kind), dimension (nx_block,ny_block,0:nvert,ngroups), intent(out) :: &
Expand Down Expand Up @@ -1762,8 +1763,8 @@ subroutine locate_triangles (nx_block, ny_block, &
ishift_br, jshift_br , & ! i,j indices of BR cell relative to edge
ishift_tc, jshift_tc , & ! i,j indices of TC cell relative to edge
ishift_bc, jshift_bc , & ! i,j indices of BC cell relative to edge
is_l, js_l , & ! i,j shifts for TL1, BL2 for area consistency
is_r, js_r , & ! i,j shifts for TR1, BR2 for area consistency
is_l, js_l , & ! i,j shifts for TL1, BL2 for area consistency
is_r, js_r , & ! i,j shifts for TR1, BR2 for area consistency
ise_tl, jse_tl , & ! i,j of TL other edge relative to edge
ise_bl, jse_bl , & ! i,j of BL other edge relative to edge
ise_tr, jse_tr , & ! i,j of TR other edge relative to edge
Expand Down Expand Up @@ -1871,7 +1872,7 @@ subroutine locate_triangles (nx_block, ny_block, &

areafac_c(:,:) = c0
areafac_ce(:,:) = c0

do ng = 1, ngroups
do j = 1, ny_block
do i = 1, nx_block
Expand Down Expand Up @@ -1908,7 +1909,7 @@ subroutine locate_triangles (nx_block, ny_block, &
ishift_bc = 0
jshift_bc = 0

! index shifts for TL1, BL2, TR1 and BR2 for area consistency
! index shifts for TL1, BL2, TR1 and BR2 for area consistency

is_l = -1
js_l = 0
Expand Down Expand Up @@ -1936,7 +1937,7 @@ subroutine locate_triangles (nx_block, ny_block, &
enddo

! area scale factor for other edge (east)

do j = 1, ny_block
do i = 1, nx_block
areafac_ce(i,j) = earea(i,j)
Expand All @@ -1960,7 +1961,7 @@ subroutine locate_triangles (nx_block, ny_block, &
ishift_bc = 0
jshift_bc = 0

! index shifts for TL1, BL2, TR1 and BR2 for area consistency
! index shifts for TL1, BL2, TR1 and BR2 for area consistency

is_l = 0
js_l = 1
Expand Down Expand Up @@ -2114,11 +2115,11 @@ subroutine locate_triangles (nx_block, ny_block, &
!-------------------------------------------------------------------
! Locate triangles in TL cell (NW for north edge, NE for east edge)
! and BL cell (W for north edge, N for east edge).
!
!
! areafact_c or areafac_ce (areafact_c for the other edge) are used
! (with shifted indices) to make sure that a flux area on one edge
! is consistent with the analogous area on the other edge and to
! ensure that areas add up when using l_fixed_area = T. See PR #849
! (with shifted indices) to make sure that a flux area on one edge
! is consistent with the analogous area on the other edge and to
! ensure that areas add up when using l_fixed_area = T. See PR #849
! for details.
!
!-------------------------------------------------------------------
Expand Down Expand Up @@ -2476,7 +2477,7 @@ subroutine locate_triangles (nx_block, ny_block, &
iflux (i,j,ng) = i + ishift_tc
jflux (i,j,ng) = j + jshift_tc
areafact(i,j,ng) = -areafac_c(i,j)

! TC2a (group 5)

ng = 5
Expand All @@ -2489,7 +2490,7 @@ subroutine locate_triangles (nx_block, ny_block, &
iflux (i,j,ng) = i + ishift_tc
jflux (i,j,ng) = j + jshift_tc
areafact(i,j,ng) = -areafac_c(i,j)

! TC3a (group 6)

ng = 6
Expand Down
25 changes: 18 additions & 7 deletions cicecore/cicedyn/general/ice_init.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -105,7 +105,7 @@ subroutine input_data
grid_ocn, grid_ocn_thrm, grid_ocn_dynu, grid_ocn_dynv, &
grid_atm, grid_atm_thrm, grid_atm_dynu, grid_atm_dynv, &
dxrect, dyrect, dxscale, dyscale, scale_dxdy, &
lonrefrect, latrefrect, pgl_global_ext
lonrefrect, latrefrect, save_ghte_ghtn
use ice_dyn_shared, only: ndte, kdyn, revised_evp, yield_curve, &
evp_algorithm, visc_method, &
seabed_stress, seabed_stress_method, &
Expand Down Expand Up @@ -375,7 +375,7 @@ subroutine input_data
ndte = 120 ! subcycles per dynamics timestep: ndte=dt_dyn/dte
evp_algorithm = 'standard_2d' ! EVP kernel (standard_2d=standard cice evp; shared_mem_1d=1d shared memory and no mpi
elasticDamp = 0.36_dbl_kind ! coefficient for calculating the parameter E
pgl_global_ext = .false. ! if true, init primary grid lengths (global ext.)
save_ghte_ghtn = .false. ! if true, save global hte and htn (global ext.)
brlx = 300.0_dbl_kind ! revised_evp values. Otherwise overwritten in ice_dyn_shared
arlx = 300.0_dbl_kind ! revised_evp values. Otherwise overwritten in ice_dyn_shared
revised_evp = .false. ! if true, use revised procedure for evp dynamics
Expand Down Expand Up @@ -962,7 +962,6 @@ subroutine input_data
call broadcast_scalar(ndte, master_task)
call broadcast_scalar(evp_algorithm, master_task)
call broadcast_scalar(elasticDamp, master_task)
call broadcast_scalar(pgl_global_ext, master_task)
call broadcast_scalar(brlx, master_task)
call broadcast_scalar(arlx, master_task)
call broadcast_scalar(revised_evp, master_task)
Expand Down Expand Up @@ -1256,6 +1255,10 @@ subroutine input_data
abort_list = trim(abort_list)//":5"
endif

if (kdyn == 1 .and. evp_algorithm == 'shared_mem_1d') then
save_ghte_ghtn = .true.
endif

if (kdyn == 2 .and. revised_evp) then
if (my_task == master_task) then
write(nu_diag,*) subname//' WARNING: revised_evp = T with EAP dynamics'
Expand Down Expand Up @@ -1294,10 +1297,10 @@ subroutine input_data
endif

if (grid_ice == 'C' .or. grid_ice == 'CD') then
if (kdyn > 1) then
if (kdyn > 1 .or. (kdyn == 1 .and. evp_algorithm /= 'standard_2d')) then
if (my_task == master_task) then
write(nu_diag,*) subname//' ERROR: grid_ice = C | CD only supported with kdyn<=1 (evp or off)'
write(nu_diag,*) subname//' ERROR: kdyn and grid_ice inconsistency'
write(nu_diag,*) subname//' ERROR: grid_ice = C | CD only supported with kdyn=1 and evp_algorithm=standard_2d'
write(nu_diag,*) subname//' ERROR: kdyn and/or evp_algorithm and grid_ice inconsistency'
endif
abort_list = trim(abort_list)//":46"
endif
Expand All @@ -1310,6 +1313,15 @@ subroutine input_data
endif
endif

if (evp_algorithm == 'shared_mem_1d' .and. &
grid_type == 'tripole') then
if (my_task == master_task) then
write(nu_diag,*) subname//' ERROR: evp_algorithm=shared_mem_1d is not tested for gridtype=tripole'
write(nu_diag,*) subname//' ERROR: change evp_algorithm to standard_2d'
endif
abort_list = trim(abort_list)//":49"
endif

capping = -9.99e30
if (kdyn == 1 .or. kdyn == 3) then
if (capping_method == 'max') then
Expand Down Expand Up @@ -1831,7 +1843,6 @@ subroutine input_data
tmpstr2 = ' : standard 2d EVP solver'
elseif (evp_algorithm == 'shared_mem_1d') then
tmpstr2 = ' : vectorized 1d EVP solver'
pgl_global_ext = .true.
else
tmpstr2 = ' : unknown value'
endif
Expand Down
2 changes: 1 addition & 1 deletion cicecore/cicedyn/general/ice_step_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -757,7 +757,7 @@ subroutine update_state (dt, daidt, dvidt, dagedt, offset)
use ice_state, only: aicen, trcrn, vicen, vsnon, &
aice, trcr, vice, vsno, aice0, trcr_depend, &
bound_state, trcr_base, nt_strata, n_trcr_strata
use ice_flux, only: Tf
use ice_flux, only: Tf
use ice_timers, only: ice_timer_start, ice_timer_stop, timer_bound, timer_updstate

real (kind=dbl_kind), intent(in) :: &
Expand Down
Loading
Loading