Skip to content

Commit

Permalink
FSD bug fixes (#424)
Browse files Browse the repository at this point in the history 
* FSD bug fixes for lateral melt and weld

* put back declare sicen

* another small welding fix

* Change print statements to use icepack_warnings

* Fix some comments

* Fix some comments

* Fix some comments

* Check for wlat present when tr_fsd = .true.

* Fix some syntax errors

---------

Co-authored-by: cmbitz <bitz@uw.edu>
  • Loading branch information
@dabail10 @cmbitz
dabail10 and cmbitz authored Mar 15, 2023
1 parent 37e215b commit a4779cc
Show file tree
Hide file tree
Showing 5 changed files with 96 additions and 79 deletions.
14 changes: 6 additions & 8 deletions columnphysics/icepack_fsd.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -699,7 +699,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 @@ -711,8 +714,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 @@ -931,7 +932,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 @@ -942,7 +942,6 @@ subroutine fsd_weld_thermo (ncat, nfsd, &
afsdn (:,:) = c0
afsd_init(:) = c0
stability = c0
prefac = p5

do n = 1, ncat

Expand Down Expand Up @@ -986,8 +985,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 @@ -1011,11 +1009,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
Original file line number Diff line number Diff line change
Expand Up @@ -884,7 +884,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 @@ -914,6 +914,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 @@ -955,7 +958,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 @@ -965,7 +968,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 @@ -979,11 +981,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 @@ -996,12 +996,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 @@ -1018,59 +1017,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 @@ -1132,8 +1129,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 @@ -1146,8 +1145,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 @@ -1243,6 +1244,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 @@ -1970,7 +1973,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 @@ -2010,10 +2013,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 @@ -2052,6 +2057,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 @@ -2122,6 +2128,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 @@ -2221,7 +2234,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

0 comments on commit a4779cc

Please sign in to comment.