Skip to content

Commit

Permalink
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge branch 'CICE-Consortium:main' into fsdfixes
Browse files Browse the repository at this point in the history
dabail10 authored Oct 25, 2024

Verified

This commit was created on GitHub.com and signed with GitHub’s verified signature. The key has expired.
2 parents a8911c3 + 6da5668 commit 07ded3b
Showing 14 changed files with 582 additions and 445 deletions.
7 changes: 5 additions & 2 deletions columnphysics/icepack_intfc.F90
Original file line number Diff line number Diff line change
@@ -112,9 +112,12 @@ module icepack_intfc
use icepack_therm_shared , only: icepack_snow_temperature
use icepack_therm_shared , only: icepack_liquidus_temperature
use icepack_therm_shared , only: icepack_sea_freezing_temperature
use icepack_therm_shared , only: icepack_init_thermo
use icepack_therm_shared , only: icepack_init_salinity
use icepack_therm_shared , only: icepack_salinity_profile
use icepack_therm_shared , only: icepack_init_trcr
use icepack_therm_shared , only: icepack_init_enthalpy
! for backwards compatibilty, remove in the future
use icepack_therm_shared , only: icepack_init_thermo => icepack_init_salinity
use icepack_therm_shared , only: icepack_init_trcr => icepack_init_enthalpy

use icepack_mushy_physics , only: icepack_enthalpy_snow
use icepack_mushy_physics , only: icepack_enthalpy_mush
39 changes: 24 additions & 15 deletions columnphysics/icepack_itd.F90
Original file line number Diff line number Diff line change
@@ -738,8 +738,8 @@ end subroutine column_conservation_check
!=======================================================================

! Cleanup subroutine that rebins thickness categories if necessary,
! eliminates very small ice areas while conserving mass and energy,
! aggregates state variables, and does a boundary call.
! eliminates very small ice areas while conserving mass and energy
! and aggregates state variables.
! It is a good idea to call this subroutine after the thermodynamics
! (thermo_vertical/thermo_itd) and again after the dynamics
! (evp/transport/ridging).
@@ -758,7 +758,8 @@ subroutine cleanup_itd (dt, hin_max, &
fpond, fresh, &
fsalt, fhocn, &
faero_ocn, fiso_ocn, &
flux_bio, Tf, limit_aice_in)
flux_bio, Tf, &
limit_aice, dorebin)

real (kind=dbl_kind), intent(in) :: &
dt ! time step
@@ -817,8 +818,9 @@ subroutine cleanup_itd (dt, hin_max, &
fiso_ocn ! isotope flux to ocean (kg/m^2/s)

logical (kind=log_kind), intent(in), optional :: &
limit_aice_in ! if false, allow aice to be out of bounds
! may want to allow this for unit tests
dorebin, & ! if false, do not call rebin (default true)
limit_aice ! if false, allow aice to be out of bounds
! may want to allow this for unit tests (default true)

! local variables

@@ -842,18 +844,25 @@ subroutine cleanup_itd (dt, hin_max, &
dflux_bio ! zapped biology flux (mmol/m^2/s)

logical (kind=log_kind) :: &
limit_aice ! if true, check for aice out of bounds
ldorebin , & ! if true, call rebin
llimit_aice ! if true, check for aice out of bounds

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

!-----------------------------------------------------------------
! Initialize
!-----------------------------------------------------------------

if (present(limit_aice_in)) then
limit_aice = limit_aice_in
if (present(limit_aice)) then
llimit_aice = limit_aice
else
limit_aice = .true.
llimit_aice = .true.
endif

if (present(dorebin)) then
ldorebin = dorebin
else
ldorebin = .true.
endif

dfpond = c0
@@ -871,7 +880,7 @@ subroutine cleanup_itd (dt, hin_max, &
call aggregate_area (aicen, aice, aice0)
if (icepack_warnings_aborted(subname)) return

if (limit_aice) then ! check for aice out of bounds
if (llimit_aice) then ! check for aice out of bounds
if (aice > c1+puny .or. aice < -puny) then
call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
call icepack_warnings_add(subname//' aggregate ice area out of bounds')
@@ -883,13 +892,13 @@ subroutine cleanup_itd (dt, hin_max, &
enddo
return
endif
endif ! limit_aice
endif ! llimit_aice

!-----------------------------------------------------------------
! Identify grid cells with ice.
!-----------------------------------------------------------------

if (aice > puny) then
if (ldorebin .and. aice > puny) then

!-----------------------------------------------------------------
! Make sure ice in each category is within its thickness bounds.
@@ -898,7 +907,7 @@ subroutine cleanup_itd (dt, hin_max, &
! correctly (e.g., very fast ice growth).
!-----------------------------------------------------------------

call rebin (trcr_depend, &
call rebin (trcr_depend, &
trcr_base, &
n_trcr_strata, &
nt_strata, &
@@ -913,7 +922,7 @@ subroutine cleanup_itd (dt, hin_max, &
! Zero out ice categories with very small areas.
!-----------------------------------------------------------------

if (limit_aice) then
if (llimit_aice) then
call zap_small_areas (dt, &
aice, aice0, &
aicen, trcrn, &
@@ -937,7 +946,7 @@ subroutine cleanup_itd (dt, hin_max, &
return
endif

endif ! l_limit_aice
endif ! llimit_aice

!-------------------------------------------------------------------
! Zap snow that has out of bounds temperatures
34 changes: 29 additions & 5 deletions columnphysics/icepack_mechred.F90
Original file line number Diff line number Diff line change
@@ -1740,7 +1740,9 @@ subroutine icepack_step_ridge(dt, ndtd, &
araftn, vraftn, &
aice, fsalt, &
first_ice, fzsal, &
flux_bio, closing, Tf )
flux_bio, closing, &
Tf, &
docleanup, dorebin)

real (kind=dbl_kind), intent(in) :: &
dt ! time step
@@ -1815,13 +1817,21 @@ subroutine icepack_step_ridge(dt, ndtd, &
logical (kind=log_kind), dimension(:), intent(inout) :: &
first_ice ! true until ice forms

logical (kind=log_kind), intent(in), optional :: &
docleanup, & ! if false, do not call cleanup_itd (default true)
dorebin ! if false, do not call rebin in cleanup_itd (default true)

!autodocument_end

! local variables

real (kind=dbl_kind) :: &
dtt ! thermo time step

logical (kind=log_kind) :: &
ldocleanup, &! if true, call cleanup_itd
ldorebin ! if true, call rebin in cleanup_itd

logical (kind=log_kind), save :: &
first_call = .true. ! first call flag

@@ -1841,6 +1851,17 @@ subroutine icepack_step_ridge(dt, ndtd, &
endif
endif

if (present(docleanup)) then
ldocleanup = docleanup
else
ldocleanup = .true.
endif

if (present(dorebin)) then
ldorebin = dorebin
else
ldorebin = .true.
endif

!-----------------------------------------------------------------
! Identify ice-ocean cells.
@@ -1880,8 +1901,9 @@ subroutine icepack_step_ridge(dt, ndtd, &
! categories with very small areas.
!-----------------------------------------------------------------

dtt = dt * ndtd ! for proper averaging over thermo timestep
call cleanup_itd (dtt, hin_max, &
if (ldocleanup) then
dtt = dt * ndtd ! for proper averaging over thermo timestep
call cleanup_itd(dtt, hin_max, &
aicen, trcrn, &
vicen, vsnon, &
aice0, aice, &
@@ -1893,8 +1915,10 @@ subroutine icepack_step_ridge(dt, ndtd, &
fpond, fresh, &
fsalt, fhocn, &
faero_ocn, fiso_ocn, &
flux_bio, Tf)
if (icepack_warnings_aborted(subname)) return
flux_bio, Tf, &
dorebin = ldorebin)
if (icepack_warnings_aborted(subname)) return
endif

first_call = .false.

8 changes: 4 additions & 4 deletions columnphysics/icepack_parameters.F90
Original file line number Diff line number Diff line change
@@ -117,8 +117,8 @@ module icepack_parameters
Timelt = 0.0_dbl_kind ,&! melting temperature, ice top surface (C)
Tsmelt = 0.0_dbl_kind ,&! melting temperature, snow top surface (C)
ice_ref_salinity =4._dbl_kind,&! (ppt)
! kice is not used for mushy thermo
kice = 2.03_dbl_kind ,&! thermal conductivity of fresh ice(W/m/deg)
! kice is only used with ktherm=1 (BL99) and conduct='MU71'
ksno = 0.30_dbl_kind ,&! thermal conductivity of snow (W/m/deg)
hs_min = 1.e-4_dbl_kind ,&! min snow thickness for computing zTsn (m)
snowpatch = 0.02_dbl_kind ,&! parameter for fractional snow area (m)
@@ -316,7 +316,7 @@ module icepack_parameters
nfreq = 25 ! number of frequencies

real (kind=dbl_kind), public :: &
floeshape = 0.66_dbl_kind ! constant from Steele (unitless)
floeshape = 0.66_dbl_kind ! constant from Rothrock 1984 (unitless)

real (kind=dbl_kind), public :: &
floediam = 300.0_dbl_kind ! effective floe diameter for lateral melt (m)
@@ -856,7 +856,7 @@ subroutine icepack_init_parameters( &
nfreq_in ! number of frequencies

real (kind=dbl_kind), intent(in), optional :: &
floeshape_in ! constant from Steele (unitless)
floeshape_in ! constant from Rothrock 1984 (unitless)

logical (kind=log_kind), intent(in), optional :: &
wave_spec_in ! if true, use wave forcing
@@ -1876,7 +1876,7 @@ subroutine icepack_query_parameters( &
nfreq_out ! number of frequencies

real (kind=dbl_kind), intent(out), optional :: &
floeshape_out ! constant from Steele (unitless)
floeshape_out ! constant from Rothrock 1984 (unitless)

logical (kind=log_kind), intent(out), optional :: &
wave_spec_out ! if true, use wave forcing
42 changes: 25 additions & 17 deletions columnphysics/icepack_therm_shared.F90
Original file line number Diff line number Diff line change
@@ -29,9 +29,9 @@ module icepack_therm_shared
public :: calculate_Tin_from_qin, &
surface_heat_flux, &
dsurface_heat_flux_dTsf, &
icepack_init_thermo, &
icepack_init_salinity, &
icepack_salinity_profile, &
icepack_init_trcr, &
icepack_init_enthalpy, &
icepack_ice_temperature, &
icepack_snow_temperature, &
icepack_liquidus_temperature, &
@@ -65,21 +65,26 @@ function calculate_Tin_from_qin (qin, Tmltk) &
Tmltk ! melting temperature at one level

real (kind=dbl_kind) :: &
Tin ! internal temperature
Tin ! internal temperature

! local variables

real (kind=dbl_kind) :: &
aa1,bb1,cc1 ! quadratic solvers
aa1,bb1,cc1,csqrt ! quadratic solvers

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

if (l_brine) then
aa1 = cp_ice
bb1 = (cp_ocn-cp_ice)*Tmltk - qin/rhoi - Lfresh
cc1 = Lfresh * Tmltk
Tin = min((-bb1 - sqrt(bb1*bb1 - c4*aa1*cc1)) / &
(c2*aa1),Tmltk)
csqrt = bb1*bb1 - c4*aa1*cc1
if (csqrt < c0) then
call icepack_warnings_add(subname//' sqrt error: ')
call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
return
endif
Tin = min((-bb1 - sqrt(csqrt)) / (c2*aa1),Tmltk)

else ! fresh ice
Tin = (Lfresh + qin/rhoi) / cp_ice
@@ -210,13 +215,14 @@ subroutine dsurface_heat_flux_dTsf(Tsf, rhoa, &
end subroutine dsurface_heat_flux_dTsf

!=======================================================================
!autodocument_start icepack_init_thermo
! Initialize the vertical profile of ice salinity and melting temperature.
!autodocument_start icepack_init_salinity
! Initialize the vertical profile of ice salinity.
! This subroutine was renamed from icepack_init_thermo in Oct 2024
!
! authors: C. M. Bitz, UW
! William H. Lipscomb, LANL

subroutine icepack_init_thermo(sprofile)
subroutine icepack_init_salinity(sprofile)

real (kind=dbl_kind), dimension(:), intent(out) :: &
sprofile ! vertical salinity profile
@@ -226,7 +232,7 @@ subroutine icepack_init_thermo(sprofile)
integer (kind=int_kind) :: k ! ice layer index
real (kind=dbl_kind) :: zn ! normalized ice thickness

character(len=*),parameter :: subname='(icepack_init_thermo)'
character(len=*),parameter :: subname='(icepack_init_salinity)'

!-----------------------------------------------------------------
! Determine l_brine based on saltmax.
@@ -239,7 +245,7 @@ subroutine icepack_init_thermo(sprofile)
if (saltmax > min_salin) l_brine = .true.

!-----------------------------------------------------------------
! Prescibe vertical profile of salinity and melting temperature.
! Prescibe vertical profile of salinity.
! Note this profile is only used for BL99 thermodynamics.
!-----------------------------------------------------------------

@@ -259,7 +265,7 @@ subroutine icepack_init_thermo(sprofile)
enddo
endif ! l_brine

end subroutine icepack_init_thermo
end subroutine icepack_init_salinity

!=======================================================================
!autodocument_start icepack_salinity_profile
@@ -282,16 +288,17 @@ function icepack_salinity_profile(zn) result(salinity)
nsal = 0.407_dbl_kind, &
msal = 0.573_dbl_kind

character(len=*),parameter :: subname='(icepack_init_thermo)'
character(len=*),parameter :: subname='(icepack_salinity_profile)'

salinity = (saltmax/c2)*(c1-cos(pi*zn**(nsal/(msal+zn))))

end function icepack_salinity_profile

!=======================================================================
!autodocument_start icepack_init_trcr
!autodocument_start icepack_init_enthalpy
! This subroutine was renamed from icepack_init_trcr in Oct 2024
!
subroutine icepack_init_trcr(Tair, Tf, &
subroutine icepack_init_enthalpy(Tair, Tf, &
Sprofile, Tprofile, &
Tsfc, &
qin, qsn)
@@ -320,7 +327,7 @@ subroutine icepack_init_trcr(Tair, Tf, &
real (kind=dbl_kind) :: &
slope, Ti

character(len=*),parameter :: subname='(icepack_init_trcr)'
character(len=*),parameter :: subname='(icepack_init_enthalpy)'

! surface temperature
Tsfc = Tf ! default
@@ -346,7 +353,7 @@ subroutine icepack_init_trcr(Tair, Tf, &
qsn(k) = -rhos*(Lfresh - cp_ice*Ti)
enddo ! nslyr

end subroutine icepack_init_trcr
end subroutine icepack_init_enthalpy

!=======================================================================
!autodocument_start icepack_liquidus_temperature
@@ -406,6 +413,7 @@ function icepack_sea_freezing_temperature(sss) result(Tf)

call icepack_warnings_add(subname//' tfrz_option unsupported: '//trim(tfrz_option))
call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
return

endif

Loading

0 comments on commit 07ded3b

Please sign in to comment.