diff --git a/src/gotm/gotm.F90 b/src/gotm/gotm.F90
index b3cf852c..a10572d2 100644
--- a/src/gotm/gotm.F90
+++ b/src/gotm/gotm.F90
@@ -784,12 +784,13 @@ subroutine integrate_gotm()
 #ifdef _ICE_
          I_0%value = transmissivity*I_0%value
          swf=melt_rate
-         shf=ocean_ice_heat_flux
-         ssf=ocean_ice_salt_flux
+         shf=ocean_ice_heat_flux + rho0*cp*T(nlev)*swf
+         ssf=ocean_ice_salt_flux + S(nlev)*swf
 #endif
       else
          swf=precip_input%value+evap
          shf=-heat_input%value ! temperature() changed to positive heat flux upwards (v7)
+         ssf= S(nlev)*swf
          tx = tx/rho0
          ty = ty/rho0
       end if
diff --git a/src/meanflow/salinity.F90 b/src/meanflow/salinity.F90
index eac8ed13..bd182cb8 100644
--- a/src/meanflow/salinity.F90
+++ b/src/meanflow/salinity.F90
@@ -116,7 +116,7 @@ subroutine salinity(nlev,dt,cnpar,wflux,sflux,nus,gams)
 !  set boundary conditions
    DiffBcup       = Neumann
    DiffBcdw       = Neumann
-   DiffSup        = -S(nlev)*wflux-sflux
+   DiffSup        = -sflux
    DiffSdw        = _ZERO_
 
    AdvBcup       = oneSided
diff --git a/src/meanflow/temperature.F90 b/src/meanflow/temperature.F90
index d314a5ad..2b767d63 100644
--- a/src/meanflow/temperature.F90
+++ b/src/meanflow/temperature.F90
@@ -137,10 +137,10 @@ subroutine temperature(nlev,dt,cnpar,I_0,wflux,hflux,nuh,gamh,rad)
    DiffBcup       = Neumann
    DiffBcdw       = Neumann
 
-!  HB:
-!  Note that this is the surface temperature flux for rigid-lid models like GOTM.
-!  For a free surface model the surface temperature flux must be - -hflux/(rho0*cp)
-   DiffTup = -T(nlev)*wflux-hflux/(rho0*cp)
+!  For the open ocean the surface temperature flux is only the diffusive
+!  component.  In case of ice cover, hflux is defined as the sum of a diffusive
+!  and advective component.
+   DiffTup = -hflux/(rho0*cp)
 #ifndef _ICE_
    ! simple sea ice model: surface heat flux switched off for sst < freezing temp
    if (T(nlev) .le. -0.0575*S(nlev)) then
diff --git a/src/util/convert_fluxes.F90 b/src/util/convert_fluxes.F90
index 6b37a80d..b2d00b05 100644
--- a/src/util/convert_fluxes.F90
+++ b/src/util/convert_fluxes.F90
@@ -66,11 +66,11 @@ subroutine  convert_fluxes(nlev,gravity,swf,shf,ssf,rad,Tsrf,Ssrf,tFlux,sFlux,bt
    beta0   = get_beta(Ssrf,Tsrf,_ZERO_)
 
    ! temperature flux and associated buoyancy flux
-   tFlux   = -Tsrf*swf - shf/(rho0*cp)
-   btFlux  =  gravity*alpha0*tFlux
+   tFlux   = - shf/(rho0*cp)
+   btFlux  = gravity*alpha0*tFlux
    
    ! salinity flux and associated buoyancy flux
-   sFlux   = -Ssrf*swf - ssf     
+   sFlux   = - ssf
    bsFlux  = -gravity*beta0*sFlux   
 
    ! radiative temperature and buoyancy flux profiles