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

FSD bug fixes #424

Merged
merged 11 commits into from
Mar 15, 2023
Merged

FSD bug fixes #424

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
5d299bd
FSD bug fixes for lateral melt and weld
cmbitz Jan 6, 2023
2cb73a0
put back declare sicen
cmbitz Jan 9, 2023
828bdc2
another small welding fix
cmbitz Jan 10, 2023
d519ed2
Change print statements to use icepack_warnings
dabail10 Jan 19, 2023
8d1e0ba
Fix some comments
dabail10 Jan 31, 2023
7be13c2
Fix some comments
dabail10 Jan 31, 2023
9632334
Fix some comments
dabail10 Jan 31, 2023
558b4b1
Fix conflict
dabail10 Jan 31, 2023
4283805
Merge branch 'CICE-Consortium:main' into FSDmods
dabail10 Feb 14, 2023
31533c7
Check for wlat present when tr_fsd = .true.
dabail10 Feb 14, 2023
e1af2cd
Fix some syntax errors
dabail10 Feb 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
7 changes: 2 additions & 5 deletions columnphysics/icepack_fsd.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -937,7 +937,6 @@ subroutine fsd_weld_thermo (ncat, nfsd, &
gain, loss ! welding tendencies

real(kind=dbl_kind) :: &
prefac , & ! multiplies kernel
kern , & ! kernel
subdt , & ! subcycling time step for stability (s)
elapsed_t ! elapsed subcycling time
Expand All @@ -948,7 +947,6 @@ subroutine fsd_weld_thermo (ncat, nfsd, &
afsdn (:,:) = c0
afsd_init(:) = c0
stability = c0
prefac = p5

do n = 1, ncat

Expand Down Expand Up @@ -992,8 +990,7 @@ subroutine fsd_weld_thermo (ncat, nfsd, &
if (k > i) then
kern = c_weld * floe_area_c(i) * aicen(n)
loss(i) = loss(i) + kern*afsd_tmp(i)*afsd_tmp(j)
if (i.eq.j) prefac = c1 ! otherwise 0.5
gain(k) = gain(k) + prefac*kern*afsd_tmp(i)*afsd_tmp(j)
gain(k) = gain(k) + kern*afsd_tmp(i)*afsd_tmp(j)
end if
end do
end do
Expand All @@ -1017,11 +1014,11 @@ subroutine fsd_weld_thermo (ncat, nfsd, &

END DO ! time

afsdn(:,n) = afsd_tmp(:)
call icepack_cleanup_fsdn (nfsd, afsdn(:,n))
if (icepack_warnings_aborted(subname)) return

do k = 1, nfsd
afsdn(k,n) = afsd_tmp(k)
trcrn(nt_fsd+k-1,n) = afsdn(k,n)
! history/diagnostics
d_afsdn_weld(k,n) = afsdn(k,n) - afsd_init(k)
Expand Down
95 changes: 43 additions & 52 deletions columnphysics/icepack_therm_itd.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -886,7 +886,7 @@ subroutine lateral_melt (dt, ncat, &
fhocn, faero_ocn, &
fiso_ocn, &
rside, meltl, &
fside, sss, &
fside, wlat, &
aicen, vicen, &
vsnon, trcrn, &
fzsal, flux_bio, &
Expand Down Expand Up @@ -916,6 +916,9 @@ subroutine lateral_melt (dt, ncat, &

real (kind=dbl_kind), intent(in) :: &
rside , & ! fraction of ice that melts laterally
wlat ! lateral melt rate (m/s)

real (kind=dbl_kind), intent(inout) :: &
fside ! lateral heat flux (W/m^2)

real (kind=dbl_kind), intent(inout) :: &
Expand Down Expand Up @@ -957,7 +960,7 @@ subroutine lateral_melt (dt, ncat, &
dfsalt , & ! change in fsalt
dvssl , & ! snow surface layer volume
dvint , & ! snow interior layer
cat1_arealoss, tmp !
bin1_arealoss, tmp !

logical (kind=log_kind) :: &
flag ! .true. if there could be lateral melting
Expand All @@ -967,7 +970,6 @@ subroutine lateral_melt (dt, ncat, &
vicen_init, & ! volume per unit area of ice (m)
G_radialn , & ! rate of lateral melt (m/s)
delta_an , & ! change in the ITD
qin , & ! enthalpy
rsiden ! delta_an/aicen

real (kind=dbl_kind), dimension (nfsd,ncat) :: &
Expand All @@ -981,11 +983,9 @@ subroutine lateral_melt (dt, ncat, &
real (kind=dbl_kind), dimension(nfsd+1) :: &
f_flx !

!echmod - for average qin
real (kind=dbl_kind), intent(in) :: &
sss
real (kind=dbl_kind) :: &
Ti, Si0, qi0, sicen, &
sicen, &
etot, & ! column energy per itd cat, for FSD code
elapsed_t, & ! FSD subcycling
subdt ! FSD timestep (s)

Expand All @@ -998,12 +998,11 @@ subroutine lateral_melt (dt, ncat, &
dfsalt = c0
dvssl = c0
dvint = c0
cat1_arealoss = c0
bin1_arealoss = c0
tmp = c0
vicen_init = c0
G_radialn = c0
delta_an = c0
qin = c0
rsiden = c0
afsdn = c1
afsdn_init = c0
Expand All @@ -1020,59 +1019,48 @@ subroutine lateral_melt (dt, ncat, &
d_afsd_latm(:) = c0
end if

if (tr_fsd .and. fside < c0) then
if (tr_fsd .and. wlat > puny) then
flag = .true.


! echmod - using category values would be preferable to the average value
! Compute average enthalpy of ice (taken from add_new_ice)
if (sss > c2 * dSin0_frazil) then
Si0 = sss - dSin0_frazil
else
Si0 = sss**2 / (c4*dSin0_frazil)
endif
Ti = min(liquidus_temperature_mush(Si0/phi_init), -p1)
qi0 = enthalpy_mush(Ti, Si0)

! for FSD rside and fside not yet computed correctly, redo here
fside = c0
do n = 1, ncat
if (ktherm == 2) then ! mushy
do k = 1, nilyr
!qin(n) = qin(n) &
! + trcrn(nt_qice+k-1,n)*vicen(n)/real(nilyr,kind=dbl_kind)
qin(n) = qi0
enddo
else
qin(n) = -rhoi*Lfresh
endif

if (qin(n) < -puny) G_radialn(n) = -fside/qin(n) ! negative
G_radialn(n) = -wlat ! negative

if (G_radialn(n) < -puny) then
if (any(afsdn(:,n) < c0)) print*,'lateral_melt B afsd < 0',n
apcraig marked this conversation as resolved.
Show resolved Hide resolved

bin1_arealoss = -trcrn(nt_fsd+1-1,n) * aicen(n) * dt &
* G_radialn(n) / floe_binwidth(1)

if (any(afsdn(:,n) < c0)) print*,&
'lateral_melt B afsd < 0',n
delta_an(n) = c0
do k = 1, nfsd
delta_an(n) = delta_an(n) + ((c2/floe_rad_c(k))*aicen(n) &
* trcrn(nt_fsd+k-1,n)*G_radialn(n)*dt) ! delta_an < 0
end do

cat1_arealoss = -trcrn(nt_fsd+1-1,n) * aicen(n) * dt &
* G_radialn(n) / floe_binwidth(1)
! add negative area loss from fsd
delta_an(n) = delta_an(n) - bin1_arealoss

delta_an(n) = c0
do k = 1, nfsd
delta_an(n) = delta_an(n) + ((c2/floe_rad_c(k))*aicen(n) &
* trcrn(nt_fsd+k-1,n)*G_radialn(n)*dt) ! delta_an < 0
end do
if (delta_an(n) > c0) print*,'ERROR delta_an > 0', delta_an(n)
apcraig marked this conversation as resolved.
Show resolved Hide resolved

! add negative area loss from fsd
delta_an(n) = delta_an(n) - cat1_arealoss
! following original code, not necessary for fsd
if (aicen(n) > c0) rsiden(n) = MIN(-delta_an(n)/aicen(n),c1)

if (delta_an(n) > c0) print*,'ERROR delta_an > 0', delta_an(n)
if (rsiden(n) < c0) print*,'ERROR rsiden < 0', rsiden(n)

! following original code, not necessary for fsd
if (aicen(n) > c0) rsiden(n) = MIN(-delta_an(n)/aicen(n),c1)
! melting energy/unit area in each column, etot < 0
etot = c0
do k = 1, nslyr
etot = etot + trcrn(nt_qsno+k-1,n) * vsnon(n)/real(nslyr,kind=dbl_kind)
enddo

if (rsiden(n) < c0) print*,'ERROR rsiden < 0', rsiden(n)
do k = 1, nilyr
etot = etot + trcrn(nt_qice+k-1,n) * vicen(n)/real(nilyr,kind=dbl_kind)
enddo ! nilyr

! lateral heat flux, fside < 0
fside = fside + rsiden(n)*etot/dt

end if ! G_radialn
enddo ! ncat

else if (rside > c0) then ! original, non-fsd implementation
Expand Down Expand Up @@ -1245,6 +1233,8 @@ subroutine lateral_melt (dt, ncat, &

if (tr_fsd) then

trcrn(nt_fsd:nt_fsd+nfsd-1,:) = afsdn

call icepack_cleanup_fsd (ncat, nfsd, trcrn(nt_fsd:nt_fsd+nfsd-1,:) )
if (icepack_warnings_aborted(subname)) return

Expand Down Expand Up @@ -1972,7 +1962,7 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, &
Tf, sss, &
salinz, &
rside, meltl, &
fside, &
fside, wlat, &
frzmlt, frazil, &
frain, fpond, &
fresh, fsalt, &
Expand Down Expand Up @@ -2012,7 +2002,7 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, &
Tf , & ! freezing temperature (C)
sss , & ! sea surface salinity (ppt)
rside , & ! fraction of ice that melts laterally
fside , & ! lateral heat flux (W/m^2)
wlat , & ! lateral melt rate
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

add units (m/s)

frzmlt , & ! freezing/melting potential (W/m^2)
wave_sig_ht ! significant height of waves in ice (m)

Expand Down Expand Up @@ -2054,6 +2044,7 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, &
real (kind=dbl_kind), intent(inout) :: &
aice , & ! sea ice concentration
aice0 , & ! concentration of open water
fside , & ! lateral heat flux (W/m^2)
frain , & ! rainfall rate (kg/m^2 s)
fpond , & ! fresh water flux to ponds (kg/m^2/s)
fresh , & ! fresh water flux to ocean (kg/m^2/s)
Expand Down Expand Up @@ -2230,7 +2221,7 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, &
fhocn, faero_ocn, &
l_fiso_ocn, &
rside, meltl, &
fside, sss, &
fside, wlat, &
aicen, vicen, &
vsnon, trcrn, &
fzsal, flux_bio, &
Expand Down
20 changes: 7 additions & 13 deletions columnphysics/icepack_therm_vertical.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -495,7 +495,7 @@ subroutine frzmlt_bottom_lateral (dt, ncat, &
strocnxT, strocnyT, &
Tbot, fbot, &
rside, Cdn_ocn, &
fside)
fside, wlat)

integer (kind=int_kind), intent(in) :: &
ncat , & ! number of thickness categories
Expand Down Expand Up @@ -529,6 +529,7 @@ subroutine frzmlt_bottom_lateral (dt, ncat, &
real (kind=dbl_kind), intent(out) :: &
Tbot , & ! ice bottom surface temperature (deg C)
fbot , & ! heat flux to ice bottom (W/m^2)
wlat , & ! lateral melt rate (m/s)
rside , & ! fraction of ice that melts laterally
fside ! lateral heat flux (W/m^2)

Expand All @@ -545,7 +546,6 @@ subroutine frzmlt_bottom_lateral (dt, ncat, &
real (kind=dbl_kind) :: &
deltaT , & ! SST - Tbot >= 0
ustar , & ! skin friction velocity for fbot (m/s)
wlat , & ! lateral melt rate (m/s)
xtmp ! temporary variable

! Parameters for bottom melting
Expand Down Expand Up @@ -618,31 +618,24 @@ subroutine frzmlt_bottom_lateral (dt, ncat, &
do n = 1, ncat

etot = c0
qavg = c0

! melting energy/unit area in each column, etot < 0

do k = 1, nslyr
etot = etot + qsnon(k,n) * vsnon(n)/real(nslyr,kind=dbl_kind)
qavg = qavg + qsnon(k,n)
enddo

do k = 1, nilyr
etot = etot + qicen(k,n) * vicen(n)/real(nilyr,kind=dbl_kind)
qavg = qavg + qicen(k,n)
enddo ! nilyr

! lateral heat flux, fside < 0
if (tr_fsd) then ! floe size distribution
fside = fside + wlat*qavg
else ! default floe size
fside = fside + rside*etot/dt
endif
fside = fside + rside*etot/dt

enddo ! n

!-----------------------------------------------------------------
! Limit bottom and lateral heat fluxes if necessary.
! FYI: fside is not yet correct for fsd, may need to adjust fbot further
!-----------------------------------------------------------------

xtmp = frzmlt/(fbot + fside + puny)
Expand Down Expand Up @@ -2121,7 +2114,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, &
fbot , &
Tbot , Tsnice , &
frzmlt , rside , &
fside , &
fside , wlat , &
fsnow , frain , &
fpond , fsloss , &
fsurf , fsurfn , &
Expand Down Expand Up @@ -2251,6 +2244,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, &
frzmlt , & ! freezing/melting potential (W/m^2)
rside , & ! fraction of ice that melts laterally
fside , & ! lateral heat flux (W/m^2)
wlat , & ! lateral melt rate (m/s)
sst , & ! sea surface temperature (C)
Tf , & ! freezing temperature (C)
Tbot , & ! ice bottom surface temperature (deg C)
Expand Down Expand Up @@ -2627,7 +2621,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, &
strocnxT, strocnyT, &
Tbot, fbot, &
rside, Cdn_ocn, &
fside)
fside, wlat)

if (icepack_warnings_aborted(subname)) return

Expand Down
1 change: 1 addition & 0 deletions configuration/driver/icedrv_flux.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -276,6 +276,7 @@ module icedrv_flux
real (kind=dbl_kind), dimension (nx), public :: &
rside , & ! fraction of ice that melts laterally
fside , & ! lateral heat flux (W/m^2)
wlat , & ! lateral heat rate (m/s)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

comment should be lateral melt rate

fsw , & ! incoming shortwave radiation (W/m^2)
coszen , & ! cosine solar zenith angle, < 0 for sun below horizon
rdg_conv, & ! convergence term for ridging (1/s)
Expand Down
6 changes: 4 additions & 2 deletions configuration/driver/icedrv_step.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -111,7 +111,7 @@ subroutine step_therm1 (dt)
use icedrv_arrays_column, only: meltsliqn, meltsliq
use icedrv_calendar, only: yday
use icedrv_domain_size, only: ncat, nilyr, nslyr, n_aero, n_iso, nx
use icedrv_flux, only: frzmlt, sst, Tf, strocnxT, strocnyT, rside, fside, &
use icedrv_flux, only: frzmlt, sst, Tf, strocnxT, strocnyT, rside, fside, wlat, &
fbot, Tbot, Tsnice
use icedrv_flux, only: meltsn, melttn, meltbn, congeln, snoicen, uatm, vatm
use icedrv_flux, only: wind, rhoa, potT, Qa, Qa_iso, zlvl, strax, stray, flatn
Expand Down Expand Up @@ -324,6 +324,7 @@ subroutine step_therm1 (dt)
fbot = fbot(i), frzmlt = frzmlt(i), &
Tbot = Tbot(i), Tsnice = Tsnice(i), &
rside = rside(i), fside = fside(i), &
wlat = wlat(i), &
fsnow = fsnow(i), frain = frain(i), &
fpond = fpond(i), fsloss = fsloss(i), &
fsurf = fsurf(i), fsurfn = fsurfn(i,:), &
Expand Down Expand Up @@ -433,7 +434,7 @@ subroutine step_therm2 (dt)
use icedrv_domain_size, only: ncat, nilyr, nslyr, n_aero, nblyr, &
nltrcr, nx, nfsd
use icedrv_flux, only: fresh, frain, fpond, frzmlt, frazil, frz_onset
use icedrv_flux, only: update_ocn_f, fsalt, Tf, sss, salinz, fhocn, rside, fside
use icedrv_flux, only: update_ocn_f, fsalt, Tf, sss, salinz, fhocn, rside, fside, wlat
use icedrv_flux, only: meltl, frazil_diag, flux_bio, faero_ocn, fiso_ocn
use icedrv_flux, only: HDO_ocn, H2_16O_ocn, H2_18O_ocn
use icedrv_init, only: tmask
Expand Down Expand Up @@ -496,6 +497,7 @@ subroutine step_therm2 (dt)
nt_strata=nt_strata(1:ntrcr,:), &
Tf=Tf(i), sss=sss(i), &
salinz=salinz(i,:), fside=fside(i), &
wlat=wlat(i), &
rside=rside(i), meltl=meltl(i), &
frzmlt=frzmlt(i), frazil=frazil(i), &
frain=frain(i), fpond=fpond(i), &
Expand Down