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

This is the code update for calling new Icepack. #215

Open
wants to merge 1 commit into
base: dev/gfdl
Choose a base branch
from
Open

This is the code update for calling new Icepack. #215

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
160 changes: 75 additions & 85 deletions config_src/external/Icepack_interfaces/icepack_mechred.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,119 +6,109 @@ module icepack_mechred
implicit none

private
public :: ridge_ice
public :: icepack_step_ridge
contains

!-----------------------------------------------------------------------
!> Interface for updating the sea-ice state due to ice ridging processes
!! using Icepack.
subroutine ridge_ice (dt, ndtd, &
ncat, n_aero, &
nilyr, nslyr, &
ntrcr, hin_max, &
rdg_conv, rdg_shear, &
aicen, trcrn, &
vicen, vsnon, &
aice0, &
trcr_depend, trcr_base, &
n_trcr_strata, &
nt_strata, &
krdg_partic, krdg_redist,&
mu_rdg, tr_brine, &
dardg1dt, dardg2dt, &
dvirdgdt, opening, &
fpond, &
fresh, fhocn, &
faero_ocn, fiso_ocn, &
aparticn, krdgn, &
aredistn, vredistn, &
dardg1ndt, dardg2ndt, &
dvirdgndt, &
araftn, vraftn, &
closing_flag,closing )

integer (kind=int_kind), intent(in) :: &
ndtd , & !< Thenumber of dynamics subcycles
ncat , & !< The number of thickness categories
nilyr , & !< The number of ice layers
nslyr , & !< The number of snow layers
n_aero, & !< The number of aerosol tracers
ntrcr !< The number of tracers in use
subroutine icepack_step_ridge(dt, ndtd, &
hin_max, &
rdg_conv, rdg_shear, &
aicen, &
trcrn, &
vicen, vsnon, &
aice0, trcr_depend, &
trcr_base, n_trcr_strata, &
nt_strata, &
dardg1dt, dardg2dt, &
dvirdgdt, opening, &
fpond, &
fresh, fhocn, &
faero_ocn, fiso_ocn, &
aparticn, krdgn, &
aredistn, vredistn, &
dardg1ndt, dardg2ndt, &
dvirdgndt, &
araftn, vraftn, &
aice, fsalt, &
first_ice, fzsal, &
flux_bio, closing, &
Tf, &
docleanup, dorebin)

real (kind=dbl_kind), intent(in) :: &
mu_rdg , & !< The e-folding scale of ridged ice [m^0.5]
dt !< The time step over which ridging occurs [s]
dt !< The time step over which ridging occurs [s]

integer (kind=int_kind), intent(in) :: &
ndtd !< The number of dynamics subcycles

real (kind=dbl_kind), dimension(0:ncat), intent(inout) :: &
hin_max !< category limits [m]

real (kind=dbl_kind), intent(in) :: &
rdg_conv , & !< Normalized energy dissipation due to convergence [s-1]
rdg_shear !< Normalized energy dissipation due to shear [s-1]

real (kind=dbl_kind), dimension (:), intent(inout) :: &
aicen , & !< The concentration of ice [nondim]
vicen , & !< The volume per unit area of ice [m]
vsnon !< The volume per unit area of snow [m]

real (kind=dbl_kind), dimension (:,:), intent(inout) :: &
trcrn !< Ice tracer concentrations [kg m-3]

real (kind=dbl_kind), intent(inout) :: &
aice0 !< The concentration of open water [nondim]

integer (kind=int_kind), dimension (:), intent(in) :: &
trcr_depend, & !< = 0 for aicen tracers, 1 for vicen, 2 for vsnon
trcr_depend, & !< = 0 for aicen tracers, 1 for vicen, 2 for vsnon
n_trcr_strata !< The number of underlying tracer layers

real (kind=dbl_kind), dimension (:,:), intent(in) :: &
trcr_base !< = 0 or 1 depending on tracer dependency [nondim]
trcr_base !< = 0 or 1 depending on tracer dependency
!! argument 2: (1) aice, (2) vice, (3) vsno [nondim]

integer (kind=int_kind), dimension (:,:), intent(in) :: &
nt_strata !< Tndices of underlying tracer layers
nt_strata !< Indices of underlying tracer layers

integer (kind=int_kind), intent(in) :: &
krdg_partic, & !< Selects participation function
krdg_redist !< Selects redistribution function

logical (kind=log_kind), intent(in) :: &
closing_flag, &!< flag if closing is valid
tr_brine !< if .true., brine height differs from ice thickness
real (kind=dbl_kind), intent(inout) :: &
aice , & !< sea ice concentration
aice0 , & !< The concentration of open water [nondim]
rdg_conv , & !< Normalized energy dissipation due to convergence [s-1]
rdg_shear, & !< Normalized energy dissipation due to shear [s-1]
dardg1dt , & !< rate of area loss by ridging ice [s-1]
dardg2dt , & !< rate of area gain by new ridges [s-1]
dvirdgdt , & !< rate of ice volume ridged [m s-1]
opening , & !< rate of opening due to divergence/shear [s-1]
fpond , & !< fresh water flux to ponds [kg m2 s-1]
fresh , & !< fresh water flux to ocean [kg m2 s-1]
fsalt , & !< salt flux to ocean [kg m2 s-1]
fhocn !< net heat flux to ocean [W m-2]

! optional history fields
real (kind=dbl_kind), intent(inout), optional :: &
dardg1dt , & !< The rate of fractional area loss by ridging ice [s-1]
dardg2dt , & !< The rate of fractional area gain by new ridges [s-1]
dvirdgdt , & !< The rate of ice volume ridged [m s-1]
opening , & !< The rate of opening due to divergence/shear [s-1]
closing , & !< The rate of closing due to divergence/shear [s-1]
fpond , & !< The fresh water flux to ponds [kg m2 s-1]
fresh , & !< The fresh water flux to ocean [kg m2 s-1]
fhocn !< The net heat flux to ocean [W m-2]
fzsal !< zsalinity flux to ocean [kg m2 s-1] (deprecated)

real (kind=dbl_kind), dimension(:), intent(inout), optional :: &
dardg1ndt , & !< The rate of fractional area loss by ridging ice [s-1]
dardg2ndt , & !< The rate of fractional area gain by new ridges [s-1]
dvirdgndt , & !< The rate of ice volume ridged [m s-1]
aparticn , & !< The participation function [nondim]
krdgn , & !< The mean ridge thickness/thickness of ridging ice [m]
araftn , & !< The rafting ice area [m2]
vraftn , & !< The rafting ice volume [m3]
aredistn , & !< The redistribution function: fraction of new ridge area [nondim]
vredistn !< The redistribution function: fraction of new ridge volume [nondim]

real (kind=dbl_kind), dimension(:), intent(inout), optional :: &
faero_ocn !< The aerosol flux to ocean [kg m-2 s-1]
real (kind=dbl_kind), intent(inout), optional :: &
closing !< rate of closing due to divergence/shear [s-1]

real (kind=dbl_kind), dimension(:), intent(inout) :: &
aicen , & !< The concentration of ice [nondim]
vicen , & !< The volume per unit area of ice [m]
vsnon , & !< The volume per unit area of snow [m]
dardg1ndt, & !< The rate of area loss by ridging ice [s-1]
dardg2ndt, & !< The rate of area gain by new ridges [s-1]
dvirdgndt, & !< The rate of ice volume ridged [m s-1]
aparticn , & !< The participation function [nondim]
krdgn , & !< The mean ridge thickness/thickness of ridging ice [m]
araftn , & !< The rafting ice area [m2]
vraftn , & !< The rafting ice volume [m3]
aredistn , & !< The redistribution function: fraction of new ridge area [nondim]
vredistn , & !< The redistribution function: fraction of new ridge volume [nondim]
faero_ocn, & !< The aerosol flux to ocean [kg m-2 s-1]
flux_bio !< All bio fluxes to ocean

real (kind=dbl_kind), dimension(:), intent(inout), optional :: &
fiso_ocn !< isotope flux to ocean [kg m-2 s-1]

fiso_ocn !< isotope flux to ocean [kg m-2 s-1]

real (kind=dbl_kind), dimension(:,:), intent(inout) :: &
trcrn !< Ice tracer concentrations [kg m-3]

logical (kind=log_kind), dimension(:), intent(inout) :: &
first_ice !< True until ice forms
real (kind=dbl_kind), intent(in) :: &
Tf !< freezing temperature
logical (kind=log_kind), dimension(:), intent(in), optional :: &
docleanup, & !< True to call cleanup_itd in Icepack
dorebin !< True to call rebin in Icepack

end subroutine ridge_ice

end subroutine icepack_step_ridge

end module icepack_mechred

Expand Down
9 changes: 9 additions & 0 deletions docs/zotero.bib
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,15 @@ @article{Hunke1997
journal = {J. Phys. Oceanography}
}

@techreport{Hunke2022,
url = {https://zenodo.org/records/7419438},
year = 2022,
title = {CICE-Consortium/Icepack: Icepack 1.3.3},
authors = {Elizabeth Hunke and Richard Allard and David A. Bailey and Philippe Blain and Anthony Craig and Frederic Dupont and
Alice DuVivier and Robert Grumbine and David Hebert and Marika Holland and Nicole Jeffery and Jean-Francois Lemieux and
Robert Osinski and Till Rasmussen and Mads Ribergaard and Lettie Roach and Andrew Roberts and Matthew Turner and Michael Winton}
}

@article{Hunke2001,
doi = {10.1006/jcph.2001.6710},
year = 2001,
Expand Down
23 changes: 21 additions & 2 deletions src/SIS_dyn_cgrid.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -160,6 +160,7 @@ module SIS_dyn_cgrid
integer :: id_sigi_hifreq = -1, id_sigii_hifreq = -1
integer :: id_stren_hifreq = -1, id_ci_hifreq = -1
integer :: id_siu = -1, id_siv = -1, id_sispeed = -1 ! SIMIP diagnostics
integer :: id_itheta = -1
!!@}
end type SIS_C_dyn_CS

Expand Down Expand Up @@ -522,6 +523,8 @@ subroutine SIS_C_dyn_init(Time, G, US, param_file, diag, CS, ntrunc)
'ice strain rate magnitude', 's-1', conversion=US%s_to_T, missing_value=missing)
CS%id_del_sh_min = register_diag_field('ice_model', 'del_sh_min', diag%axesT1, Time, &
'minimum ice strain rate magnitude', 's-1', conversion=US%s_to_T, missing_value=missing)
CS%id_itheta = register_diag_field('ice_model', 'itheta', diag%axesT1, Time, &
'ice atan(shear/divergence)', 'rad', missing_value=missing)

CS%id_ui_hifreq = register_diag_field('ice_model', 'ui_hf', diag%axesCu1, Time, &
'ice velocity - x component', 'm/s', missing_value=missing, &
Expand Down Expand Up @@ -637,6 +640,7 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, fxat, fyat, &
dx2T, dy2T, & ! dx^2 or dy^2 at T points [L2 ~> m2].
dx_dyT, dy_dxT, & ! dx/dy or dy_dx at T points [nondim].
siu, siv, sispeed, & ! diagnostics on T points [L T-1 ~> m s-1].
itheta, & ! Angle given by atan(shear/divergence)
ui_east, & ! Surface velocity due east component [L T-1 ~> m s-1]
vi_north ! Surface velocity due north component [L T-1 ~> m s-1]

Expand Down Expand Up @@ -745,6 +749,7 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, fxat, fyat, &
real :: m_neglect2 ! A tiny mass per unit area squared [R2 Z2 ~> kg2 m-4].
real :: m_neglect4 ! A tiny mass per unit area to the 4th power [R4 Z4 ~> kg4 m-8].
real :: sum_area ! The sum of ocean areas around a vorticity point [L2 ~> m2].
real :: half_pi ! pi/2.

type(time_type) :: &
time_it_start, & ! The starting time of the iterative steps.
Expand All @@ -764,6 +769,7 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, fxat, fyat, &
"SIS_C_dynamics is written to require a 2-point halo or 1-point and symmetric memory.")

halo_sh_Ds = min(isc-G%isd, jsc-G%jsd, 2)
half_pi = 2 * atan(1.0)

! Zero these arrays to accumulate sums.
fxoc(:,:) = 0.0 ; fyoc(:,:) = 0.0
Expand Down Expand Up @@ -1062,14 +1068,26 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, fxat, fyat, &

! calculate viscosities - how often should we do this ?
!$OMP parallel do default(none) shared(isc,iec,jsc,jec,del_sh,zeta,sh_Dd,sh_Dt, &
!$OMP I_EC2,sh_Ds,pres_mice,mice,del_sh_min_pr)
!$OMP I_EC2,sh_Ds,pres_mice,mice,del_sh_min_pr, &
!$OMP itheta)
if (CS%id_itheta > 0) then
do j=jsc-1,jec+1 ; do i=isc-1,iec+1
itheta(i,j) = 0.0
enddo ; enddo
endif
do j=jsc-1,jec+1 ; do i=isc-1,iec+1
! Averaging the squared shearing strain is larger than squaring
! the averaged strain. I don't know what is better. -RWH
del_sh(i,j) = sqrt(sh_Dd(i,j)**2 + I_EC2 * (sh_Dt(i,j)**2 + &
(0.25 * ((sh_Ds(I-1,J-1) + sh_Ds(I,J)) + &
(sh_Ds(I-1,J) + sh_Ds(I,J-1))))**2 ) ) ! H&D eqn 9

if (CS%id_itheta > 0 .and. ci(i,j) > 0.0 .and. sh_Dd(i,j) /= 0.0) then
itheta(i,j) = atan( 0.25 * ((sh_Ds(I-1,J-1) + sh_Ds(I,J)) + &
(sh_Ds(I-1,J) + sh_Ds(I,J-1))) &
/ abs(sh_Dd(i,j)) )
if (itheta(i,j) < 0.0) itheta(i,j) = itheta(i,j) + half_pi
endif
if (max(del_sh(i,j), del_sh_min_pr(i,j)*pres_mice(i,j)) /= 0.) then
zeta(i,j) = 0.5*pres_mice(i,j)*mice(i,j) / &
max(del_sh(i,j), del_sh_min_pr(i,j)*pres_mice(i,j))
Expand Down Expand Up @@ -1555,6 +1573,7 @@ subroutine SIS_C_dynamics(ci, mis, mice, ui, vi, uo, vo, fxat, fyat, &
if (CS%id_sh_s>0) call post_SIS_data(CS%id_sh_s, sh_Ds, CS%diag)

if (CS%id_del_sh>0) call post_SIS_data(CS%id_del_sh, del_sh, CS%diag)
if (CS%id_itheta>0) call post_SIS_data(CS%id_itheta, itheta, CS%diag)
if (CS%id_del_sh_min>0) then
do j=jsc,jec ; do i=isc,iec
diag_val(i,j) = del_sh_min_pr(i,j)*pres_mice(i,j)
Expand Down Expand Up @@ -1652,7 +1671,7 @@ subroutine limit_stresses(pres_mice, mice, str_d, str_t, str_s, G, US, CS, limit
enddo ; enddo

! This commented out version seems to work, but is not obviously better than
! treating each component separately, and the later is simpler.
! treating each component separately, and the latter is simpler.
! EC2 = CS%EC**2
! do J=jsc-1,jec ; do I=isc-1,iec
! ! Rescale str_s based on interpolated values of str_d and str_t, which works
Expand Down
18 changes: 10 additions & 8 deletions src/SIS_dyn_trans.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -407,7 +407,7 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, U
G, US, IG, CS, OBC)

! Complete the category-resolved mass and tracer transport and update the ice state type.
call complete_IST_transport(CS%DS2d, CS%CAS, IST, dt_adv_cycle, G, US, IG, CS)
call complete_IST_transport(CS%DS2d, CS%CAS, IST, OSS, dt_adv_cycle, G, US, IG, CS)

else ! (.not.CS%merged_cont)

Expand Down Expand Up @@ -624,10 +624,10 @@ subroutine SIS_dynamics_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G, U
if (CS%do_ridging) then
call finish_ice_transport(CS%CAS, IST, IST%TrReg, G, US, IG, dt_slow, CS%SIS_transport_CSp, &
! rdg_rate=DS2d%avg_ridge_rate)
rdg_rate=IST%rdg_rate)
OSS=OSS, rdg_rate=IST%rdg_rate)
DS2d%ridge_rate_count = 0. ; DS2d%avg_ridge_rate(:,:) = 0.0
else
call finish_ice_transport(CS%CAS, IST, IST%TrReg, G, US, IG, dt_slow,CS%SIS_transport_CSp)
call finish_ice_transport(CS%CAS, IST, IST%TrReg, G, US, IG, dt_slow, CS%SIS_transport_CSp)
endif
endif
call cpu_clock_end(iceClock8)
Expand Down Expand Up @@ -716,7 +716,7 @@ subroutine SIS_multi_dyn_trans(IST, OSS, FIA, IOF, dt_slow, CS, icebergs_CS, G,
! Complete the category-resolved mass and tracer transport and update the ice state type.
! This must be done before the next thermodynamic step.
if (end_of_cycle) &
call complete_IST_transport(CS%DS2d, CS%CAS, IST, dt_adv_cycle, G, US, IG, CS)
call complete_IST_transport(CS%DS2d, CS%CAS, IST, OSS, dt_adv_cycle, G, US, IG, CS)

if (CS%column_check .and. IST%valid_IST) & ! This is just here from early debugging exercises,
call write_ice_statistics(IST, CS%Time, CS%n_calls, G, US, IG, CS%sum_output_CSp, &
Expand All @@ -735,14 +735,16 @@ end subroutine SIS_multi_dyn_trans

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
!> Complete the category-resolved mass and tracer transport and update the ice state type.
subroutine complete_IST_transport(DS2d, CAS, IST, dt_adv_cycle, G, US, IG, CS)
subroutine complete_IST_transport(DS2d, CAS, IST, OSS, dt_adv_cycle, G, US, IG, CS)
type(ice_state_type), intent(inout) :: IST !< A type describing the state of the sea ice
type(ocean_sfc_state_type), intent(in) :: OSS !< A structure containing the arrays that describe
!! the ocean's surface state for the ice model.
type(dyn_state_2d), intent(inout) :: DS2d !< A simplified 2-d description of the ice state
!! integrated across thickness categories and layers.
!! integrated across thickness categories and layers.
type(cell_average_state_type), intent(inout) :: CAS !< A structure with ocean-cell averaged masses.
real, intent(in) :: dt_adv_cycle !< The time since the last IST transport [T ~> s].
type(SIS_hor_grid_type), intent(inout) :: G !< The horizontal grid type
type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors
type(unit_scale_type), intent(in) :: US !< A structure with unit conversion factors
type(ice_grid_type), intent(inout) :: IG !< The sea-ice specific grid type
type(dyn_trans_CS), pointer :: CS !< The control structure for the SIS_dyn_trans module

Expand All @@ -763,7 +765,7 @@ subroutine complete_IST_transport(DS2d, CAS, IST, dt_adv_cycle, G, US, IG, CS)
if (CS%do_ridging) then
call finish_ice_transport(CS%CAS, IST, IST%TrReg, G, US, IG, dt_adv_cycle, CS%SIS_transport_CSp, &
! rdg_rate=DS2d%avg_ridge_rate)
rdg_rate=IST%rdg_rate)
OSS=OSS, rdg_rate=IST%rdg_rate)
DS2d%ridge_rate_count = 0. ; DS2d%avg_ridge_rate(:,:) = 0.0
else
call finish_ice_transport(CS%CAS, IST, IST%TrReg, G, US, IG, dt_adv_cycle, CS%SIS_transport_CSp)
Expand Down
Loading