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

Update to CICE-Consortium main, CICE6.4.1 release, #f813294, Dec 9, 2022 #106

Merged
merged 5 commits into from
Dec 16, 2022
Merged

Update to CICE-Consortium main, CICE6.4.1 release, #f813294, Dec 9, 2022 #106

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
b16d7fd
Update Icepack to #5999551 including snowbrinebugs (#797)
apcraig Dec 2, 2022
4befa1e
Saltflux option for CICE (#799)
dabail10 Dec 7, 2022
48cf07a
Update stress12T calculation for C-grid (#802)
apcraig Dec 7, 2022
8df2d7b
Merge branch 'main' of https://github.com/apcraig/cice into e3sm_upd641
apcraig Dec 15, 2022
b28511b
Update to E3SM-Project/Icepack/cice-consortium/E3SM-icepack-initial-i…
apcraig Dec 15, 2022
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
45 changes: 37 additions & 8 deletions cicecore/cicedyn/analysis/ice_diagnostics.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -87,6 +87,8 @@ module ice_diagnostics
totms , & ! total ice/snow water mass (sh)
totmin , & ! total ice water mass (nh)
totmis , & ! total ice water mass (sh)
totsn , & ! total salt mass (nh)
totss , & ! total salt mass (sh)
toten , & ! total ice/snow energy (J)
totes ! total ice/snow energy (J)

Expand Down Expand Up @@ -154,14 +156,16 @@ subroutine runtime_diags (dt)
rhofresh, lfresh, lvap, ice_ref_salinity, Tffresh

character (len=char_len) :: &
snwredist
snwredist, saltflux_option

! hemispheric state quantities
real (kind=dbl_kind) :: &
umaxn, hmaxn, shmaxn, arean, snwmxn, extentn, shmaxnt, &
umaxs, hmaxs, shmaxs, areas, snwmxs, extents, shmaxst, &
etotn, mtotn, micen, msnwn, pmaxn, ketotn, &
etots, mtots, mices, msnws, pmaxs, ketots, &
stotn, &
stots, &
urmsn, albtotn, arean_alb, mpndn, ptotn, spondn, &
urmss, albtots, areas_alb, mpnds, ptots, sponds

Expand Down Expand Up @@ -226,7 +230,7 @@ subroutine runtime_diags (dt)
awtvdr_out=awtvdr, awtidr_out=awtidr, awtvdf_out=awtvdf, awtidf_out=awtidf, &
rhofresh_out=rhofresh, lfresh_out=lfresh, lvap_out=lvap, &
ice_ref_salinity_out=ice_ref_salinity,snwredist_out=snwredist, &
snwgrain_out=snwgrain)
snwgrain_out=snwgrain, saltflux_option_out=saltflux_option)
call icepack_warnings_flush(nu_diag)
if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
file=__FILE__, line=__LINE__)
Expand Down Expand Up @@ -512,6 +516,15 @@ subroutine runtime_diags (dt)
etots = global_sum(work1, distrb_info, &
field_loc_center, tareas)

! total salt volume
call total_salt (work2)

stotn = global_sum(work2, distrb_info, &
field_loc_center, tarean)
stots = global_sum(work2, distrb_info, &
field_loc_center, tareas)


!-----------------------------------------------------------------
! various fluxes
!-----------------------------------------------------------------
Expand Down Expand Up @@ -785,12 +798,22 @@ subroutine runtime_diags (dt)
swerrs = (fswnets - fswdns) / (fswnets - c1)

! salt mass
msltn = micen*ice_ref_salinity*p001
mslts = mices*ice_ref_salinity*p001
if (saltflux_option == 'prognostic') then
! compute the total salt mass
msltn = stotn*rhoi*p001
mslts = stots*rhoi*p001

! change in salt mass
delmsltn = rhoi*(stotn-totsn)*p001
delmslts = rhoi*(stots-totss)*p001
else
msltn = micen*ice_ref_salinity*p001
mslts = mices*ice_ref_salinity*p001

! change in salt mass
delmsltn = delmxn*ice_ref_salinity*p001
delmslts = delmxs*ice_ref_salinity*p001
! change in salt mass
delmsltn = delmxn*ice_ref_salinity*p001
delmslts = delmxs*ice_ref_salinity*p001
endif

! salt error
serrn = (sfsaltn + delmsltn) / (msltn + c1)
Expand Down Expand Up @@ -1275,7 +1298,7 @@ subroutine init_mass_diags
rhoi, rhos, rhofresh

real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
work1
work1, work2

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

Expand Down Expand Up @@ -1310,6 +1333,12 @@ subroutine init_mass_diags
toten = global_sum(work1, distrb_info, field_loc_center, tarean)
totes = global_sum(work1, distrb_info, field_loc_center, tareas)

! north/south salt
call total_salt (work2)
totsn = global_sum(work2, distrb_info, field_loc_center, tarean)
totss = global_sum(work2, distrb_info, field_loc_center, tareas)


if (print_points) then
do n = 1, npnt
if (my_task == pmloc(n)) then
Expand Down
21 changes: 14 additions & 7 deletions cicecore/cicedyn/analysis/ice_history.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -420,10 +420,6 @@ subroutine init_hist (dt)
f_taubyE = f_tauby
endif

#ifndef ncdf
f_bounds = .false.
#endif

! write dimensions for 3D or 4D history variables
! note: list of variables checked here is incomplete
if (f_aicen(1:1) /= 'x' .or. f_vicen(1:1) /= 'x' .or. &
Expand Down Expand Up @@ -2168,12 +2164,13 @@ subroutine accum_hist (dt)

real (kind=dbl_kind) :: awtvdr, awtidr, awtvdf, awtidf, puny, secday, rad_to_deg
real (kind=dbl_kind) :: Tffresh, rhoi, rhos, rhow, ice_ref_salinity
real (kind=dbl_kind) :: rho_ice, rho_ocn, Tice, Sbr, phi, rhob, dfresh, dfsalt
real (kind=dbl_kind) :: rho_ice, rho_ocn, Tice, Sbr, phi, rhob, dfresh, dfsalt, sicen
logical (kind=log_kind) :: formdrag, skl_bgc
logical (kind=log_kind) :: tr_pond, tr_aero, tr_brine, tr_snow
integer (kind=int_kind) :: ktherm
integer (kind=int_kind) :: nt_sice, nt_qice, nt_qsno, nt_iage, nt_FY, nt_Tsfc, &
nt_alvl, nt_vlvl
character (len=char_len) :: saltflux_option

type (block) :: &
this_block ! block information for current block
Expand All @@ -2185,6 +2182,7 @@ subroutine accum_hist (dt)
call icepack_query_parameters(Tffresh_out=Tffresh, rhoi_out=rhoi, rhos_out=rhos, &
rhow_out=rhow, ice_ref_salinity_out=ice_ref_salinity)
call icepack_query_parameters(formdrag_out=formdrag, skl_bgc_out=skl_bgc, ktherm_out=ktherm)
call icepack_query_parameters(saltflux_option_out=saltflux_option)
call icepack_query_tracer_flags(tr_pond_out=tr_pond, tr_aero_out=tr_aero, &
tr_brine_out=tr_brine, tr_snow_out=tr_snow)
call icepack_query_tracer_indices(nt_sice_out=nt_sice, nt_qice_out=nt_qice, &
Expand Down Expand Up @@ -2269,7 +2267,7 @@ subroutine accum_hist (dt)
!---------------------------------------------------------------

!$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block, &
!$OMP k,n,qn,ns,sn,rho_ocn,rho_ice,Tice,Sbr,phi,rhob,dfresh,dfsalt, &
!$OMP k,n,qn,ns,sn,rho_ocn,rho_ice,Tice,Sbr,phi,rhob,dfresh,dfsalt,sicen, &
!$OMP worka,workb,worka3,Tinz4d,Sinz4d,Tsnz4d)

do iblk = 1, nblocks
Expand Down Expand Up @@ -3228,7 +3226,16 @@ subroutine accum_hist (dt)
dfresh = -rhoi*frazil(i,j,iblk)/dt
endif
endif
dfsalt = ice_ref_salinity*p001*dfresh
if (saltflux_option == 'prognostic') then
sicen = c0
do k = 1, nzilyr
sicen = sicen + trcr(i,j,nt_sice+k-1,iblk)*vice(i,j,iblk) &
/ real(nzilyr,kind=dbl_kind)
enddo
dfsalt = sicen*p001*dfresh
else
dfsalt = ice_ref_salinity*p001*dfresh
endif
worka(i,j) = aice(i,j,iblk)*(fsalt(i,j,iblk)+dfsalt)
endif
enddo
Expand Down
2 changes: 1 addition & 1 deletion cicecore/cicedyn/analysis/ice_history_shared.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@ module ice_history_shared
logical (kind=log_kind), public :: &
hist_avg ! if true, write averaged data instead of snapshots

character (len=char_len), public :: &
character (len=char_len_long), public :: &
history_file , & ! output file for history
incond_file ! output file for snapshot initial conditions

Expand Down
89 changes: 52 additions & 37 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 @@ -945,7 +945,8 @@ subroutine evp (dt)
uarea (:,:,iblk), DminTarea (:,:,iblk), &
strength (:,:,iblk), shearU (:,:,iblk), &
zetax2T (:,:,iblk), etax2T (:,:,iblk), &
stresspT (:,:,iblk), stressmT (:,:,iblk))
stresspT (:,:,iblk), stressmT (:,:,iblk), &
stress12T (:,:,iblk))

enddo
!$OMP END PARALLEL DO
Expand Down Expand Up @@ -1730,7 +1731,8 @@ subroutine stressC_T (nx_block, ny_block , &
uarea , DminTarea , &
strength , shearU , &
zetax2T , etax2T , &
stressp , stressm )
stresspT , stressmT , &
stress12T)

use ice_dyn_shared, only: strain_rates_T, capping, &
visc_replpress, e_factor
Expand All @@ -1744,37 +1746,40 @@ subroutine stressC_T (nx_block, ny_block , &
indxTj ! compressed index in j-direction

real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
uvelE , & ! x-component of velocity (m/s) at the E point
vvelE , & ! y-component of velocity (m/s) at the E point
uvelN , & ! x-component of velocity (m/s) at the N point
vvelN , & ! y-component of velocity (m/s) at the N point
dxN , & ! width of N-cell through the middle (m)
dyE , & ! height of E-cell through the middle (m)
dxT , & ! width of T-cell through the middle (m)
dyT , & ! height of T-cell through the middle (m)
strength , & ! ice strength (N/m)
shearU , & ! shearU local for this routine
uarea , & ! area of u cell
DminTarea ! deltaminEVP*tarea
uvelE , & ! x-component of velocity (m/s) at the E point
vvelE , & ! y-component of velocity (m/s) at the E point
uvelN , & ! x-component of velocity (m/s) at the N point
vvelN , & ! y-component of velocity (m/s) at the N point
dxN , & ! width of N-cell through the middle (m)
dyE , & ! height of E-cell through the middle (m)
dxT , & ! width of T-cell through the middle (m)
dyT , & ! height of T-cell through the middle (m)
strength , & ! ice strength (N/m)
shearU , & ! shearU local for this routine
uarea , & ! area of u cell
DminTarea ! deltaminEVP*tarea

real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
zetax2T , & ! zetax2 = 2*zeta (bulk viscosity)
etax2T , & ! etax2 = 2*eta (shear viscosity)
stressp , & ! sigma11+sigma22
stressm ! sigma11-sigma22
zetax2T , & ! zetax2 = 2*zeta (bulk viscosity)
etax2T , & ! etax2 = 2*eta (shear viscosity)
stresspT , & ! sigma11+sigma22
stressmT , & ! sigma11-sigma22
stress12T ! sigma12

! local variables

integer (kind=int_kind) :: &
i, j, ij

real (kind=dbl_kind), dimension (nx_block,ny_block) :: &
divT , & ! divergence at T point
tensionT ! tension at T point
divT , & ! divergence at T point
tensionT ! tension at T point

real (kind=dbl_kind) :: &
shearTsqr , & ! strain rates squared at T point
shearT , & ! strain rate at T point
DeltaT , & ! delt at T point
uareaavgr , & ! 1 / uarea avg
rep_prsT ! replacement pressure at T point

character(len=*), parameter :: subname = '(stressC_T)'
Expand All @@ -1801,11 +1806,19 @@ subroutine stressC_T (nx_block, ny_block , &
! U point values (Bouillon et al., 2013, Kimmritz et al., 2016
!-----------------------------------------------------------------

uareaavgr = c1/(uarea(i,j)+uarea(i,j-1)+uarea(i-1,j-1)+uarea(i-1,j))

shearTsqr = (shearU(i ,j )**2 * uarea(i ,j ) &
+ shearU(i ,j-1)**2 * uarea(i ,j-1) &
+ shearU(i-1,j-1)**2 * uarea(i-1,j-1) &
+ shearU(i-1,j )**2 * uarea(i-1,j )) &
/ (uarea(i,j)+uarea(i,j-1)+uarea(i-1,j-1)+uarea(i-1,j))
* uareaavgr

shearT = (shearU(i ,j ) * uarea(i ,j ) &
+ shearU(i ,j-1) * uarea(i ,j-1) &
+ shearU(i-1,j-1) * uarea(i-1,j-1) &
+ shearU(i-1,j ) * uarea(i-1,j )) &
* uareaavgr

DeltaT = sqrt(divT(i,j)**2 + e_factor*(tensionT(i,j)**2 + shearTsqr))

Expand All @@ -1822,11 +1835,14 @@ subroutine stressC_T (nx_block, ny_block , &

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

stressp(i,j) = (stressp(i,j)*(c1-arlx1i*revp) &
+ arlx1i*(zetax2T(i,j)*divT(i,j) - rep_prsT)) * denom1
stresspT(i,j) = (stresspT (i,j)*(c1-arlx1i*revp) &
+ arlx1i*(zetax2T(i,j)*divT(i,j) - rep_prsT)) * denom1

stressm(i,j) = (stressm(i,j)*(c1-arlx1i*revp) &
+ arlx1i*etax2T(i,j)*tensionT(i,j)) * denom1
stressmT(i,j) = (stressmT (i,j)*(c1-arlx1i*revp) &
+ arlx1i*etax2T(i,j)*tensionT(i,j)) * denom1

stress12T(i,j) = (stress12T(i,j)*(c1-arlx1i*revp) &
+ arlx1i*p5*etax2T(i,j)*shearT ) * denom1

enddo ! ij

Expand All @@ -1851,7 +1867,7 @@ subroutine stressC_U (nx_block , ny_block ,&
uarea , &
etax2U , deltaU ,&
strengthU, shearU ,&
stress12)
stress12U)

use ice_dyn_shared, only: visc_replpress, &
visc_method, deltaminEVP, capping
Expand All @@ -1872,7 +1888,7 @@ subroutine stressC_U (nx_block , ny_block ,&
strengthU ! ice strength at the U point

real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
stress12 ! sigma12
stress12U ! sigma12

! local variables

Expand All @@ -1891,15 +1907,15 @@ subroutine stressC_U (nx_block , ny_block ,&
! viscosities and replacement pressure at U point
! avg_zeta: Bouillon et al. 2013, C1 method of Kimmritz et al. 2016
! avg_strength: C2 method of Kimmritz et al. 2016
! if outside do and stress12 equation repeated in each loop for performance
! if outside do and stress12U equation repeated in each loop for performance
!-----------------------------------------------------------------

if (visc_method == 'avg_zeta') then
do ij = 1, icellU
i = indxUi(ij)
j = indxUj(ij)
stress12(i,j) = (stress12(i,j)*(c1-arlx1i*revp) &
+ arlx1i*p5*etax2U(i,j)*shearU(i,j)) * denom1
stress12U(i,j) = (stress12U(i,j)*(c1-arlx1i*revp) &
+ arlx1i*p5*etax2U(i,j)*shearU(i,j)) * denom1
enddo

elseif (visc_method == 'avg_strength') then
Expand All @@ -1911,8 +1927,8 @@ subroutine stressC_U (nx_block , ny_block ,&
! minimal extra calculations here even though it seems like there is
call visc_replpress (strengthU(i,j), DminUarea, deltaU(i,j), &
lzetax2U , letax2U , lrep_prsU , capping)
stress12(i,j) = (stress12(i,j)*(c1-arlx1i*revp) &
+ arlx1i*p5*letax2U*shearU(i,j)) * denom1
stress12U(i,j) = (stress12U(i,j)*(c1-arlx1i*revp) &
+ arlx1i*p5*letax2U*shearU(i,j)) * denom1
enddo

endif
Expand Down Expand Up @@ -1976,7 +1992,7 @@ subroutine stressCD_T (nx_block, ny_block , &
real (kind=dbl_kind), dimension (nx_block,ny_block) :: &
divT , & ! divergence at T point
tensionT , & ! tension at T point
shearT , & ! sheat at T point
shearT , & ! shear at T point
DeltaT ! delt at T point

real (kind=dbl_kind) :: &
Expand All @@ -1991,7 +2007,7 @@ subroutine stressCD_T (nx_block, ny_block , &

call strain_rates_T (nx_block , ny_block , &
icellT , &
indxTi (:), indxTj (:) , &
indxTi(:) , indxTj (:) , &
uvelE (:,:), vvelE (:,:), &
uvelN (:,:), vvelN (:,:), &
dxN (:,:), dyE (:,:), &
Expand Down Expand Up @@ -2046,8 +2062,7 @@ subroutine stressCD_U (nx_block, ny_block, &
stresspU , stressmU, &
stress12U)

use ice_dyn_shared, only: strain_rates_U, &
visc_replpress, &
use ice_dyn_shared, only: visc_replpress, &
visc_method, deltaminEVP, capping

integer (kind=int_kind), intent(in) :: &
Expand All @@ -2059,7 +2074,7 @@ subroutine stressCD_U (nx_block, ny_block, &
indxUj ! compressed index in j-direction

real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
uarea , & ! area of U-cell (m^2)
uarea , & ! area of U-cell (m^2)
zetax2U , & ! 2*zeta at U point
etax2U , & ! 2*eta at U point
strengthU, & ! ice strength at U point
Expand Down
2 changes: 1 addition & 1 deletion cicecore/cicedyn/dynamics/ice_dyn_evp_1d.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -889,7 +889,7 @@ subroutine evp1d_halo_update(NAVEL_len, lb, ub, uvel, vvel, &

#ifdef _OPENACC
!$acc parallel &
!$acc present(uvel, vvel) &
!$acc present(uvel, vvel)
!$acc loop
do iw = 1, NAVEL_len
if (halo_parent(iw) == 0) cycle
Expand Down
Loading