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 all 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
14 changes: 6 additions & 8 deletions columnphysics/icepack_fsd.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -705,7 +705,10 @@ subroutine fsd_add_new_ice (ncat, n, nfsd, &
DO WHILE (elapsed_t.lt.dt)

nsubt = nsubt + 1
if (nsubt.gt.100) print *, 'latg not converging'
if (nsubt.gt.100) then
write(warnstr,*) subname,'latg not converging'
call icepack_warnings_add(warnstr)
endif

! finite differences
df_flx(:) = c0 ! NB could stay zero if all in largest FS cat
Expand All @@ -717,8 +720,6 @@ subroutine fsd_add_new_ice (ncat, n, nfsd, &
df_flx(k) = f_flx(k+1) - f_flx(k)
end do

! if (abs(sum(df_flx)) > puny) print*,'fsd_add_new ERROR df_flx /= 0'

dafsd_tmp(:) = c0
do k = 1, nfsd
dafsd_tmp(k) = (-df_flx(k) + c2 * G_radial * afsdn_latg(k,n) &
Expand Down Expand Up @@ -937,7 +938,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 +948,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 +991,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 +1015,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
125 changes: 69 additions & 56 deletions columnphysics/icepack_therm_itd.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -887,7 +887,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 @@ -917,6 +917,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 @@ -958,7 +961,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 @@ -968,7 +971,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 @@ -982,11 +984,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 @@ -999,12 +999,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 @@ -1021,59 +1020,57 @@ 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)) then
write(warnstr,*) subname, 'lateral_melt B afsd < 0 ',n
call icepack_warnings_add(warnstr)
endif

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) then
write(warnstr,*) subname, 'ERROR delta_an > 0 ',delta_an(n)
call icepack_warnings_add(warnstr)
endif

! 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) then
write(warnstr,*) subname, 'ERROR rsiden < 0 ',rsiden(n)
call icepack_warnings_add(warnstr)
endif

! 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 @@ -1135,8 +1132,10 @@ subroutine lateral_melt (dt, ncat, &
DO WHILE (elapsed_t.lt.dt)

nsubt = nsubt + 1
if (nsubt.gt.100) &
print *, 'latm not converging'
if (nsubt.gt.100) then
write(warnstr,*) subname, 'latm not converging'
call icepack_warnings_add(warnstr)
endif

! finite differences
df_flx(:) = c0
Expand All @@ -1149,8 +1148,10 @@ subroutine lateral_melt (dt, ncat, &
df_flx(k) = f_flx(k+1) - f_flx(k)
end do

if (abs(sum(df_flx(:))) > puny) &
print*,'sum(df_flx)/=0'
if (abs(sum(df_flx(:))) > puny) then
write(warnstr,*) subname, 'sum(df_flx) /= 0'
call icepack_warnings_add(warnstr)
endif

! this term ensures area conservation
tmp = SUM(afsd_tmp(:)/floe_rad_c(:))
Expand Down Expand Up @@ -1246,6 +1247,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 @@ -1973,7 +1976,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 @@ -2013,10 +2016,12 @@ 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)
frzmlt , & ! freezing/melting potential (W/m^2)
wave_sig_ht ! significant height of waves in ice (m)

real (kind=dbl_kind), intent(in), optional :: &
wlat ! lateral melt rate (m/s)

real (kind=dbl_kind), dimension(:), intent(in) :: &
wave_spectrum ! ocean surface wave spectrum E(f) (m^2 s)

Expand Down Expand Up @@ -2055,6 +2060,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 @@ -2125,6 +2131,13 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, &
return
endif
endif
if (tr_fsd) then
if (.not.(present(wlat))) then
call icepack_warnings_add(subname//' error in FSD arguments, tr_fsd=T')
call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
return
endif
endif
endif

!-----------------------------------------------------------------
Expand Down Expand Up @@ -2224,7 +2237,7 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr, &
fhocn, faero_ocn, &
fiso_ocn, &
rside, meltl, &
fside, sss, &
fside, wlat, &
aicen, vicen, &
vsnon, trcrn, &
fzsal, flux_bio, &
Expand Down
Loading