From d5cdca1fe442c9e3ea888bda79854357784cf278 Mon Sep 17 00:00:00 2001
From: Ou Wang <owang01@yahoo.com>
Date: Wed, 3 Apr 2024 17:06:59 -0700
Subject: [PATCH] Implement sea-ice freeboard cost

---
 pkg/ecco/cost_gencost_customize.F |  3 +++
 pkg/ecco/ecco_phys.F              | 41 +++++++++++++++++++++++++++++--
 2 files changed, 42 insertions(+), 2 deletions(-)

diff --git a/pkg/ecco/cost_gencost_customize.F b/pkg/ecco/cost_gencost_customize.F
index 303b135922..24598e42bb 100644
--- a/pkg/ecco/cost_gencost_customize.F
+++ b/pkg/ecco/cost_gencost_customize.F
@@ -221,6 +221,9 @@ subroutine cost_gencost_customize( mythid )
          elseif (gencost_barfile(k)(1:9).EQ.'m_sihsnow') then
            gencost_modfld(i,j,bi,bj,k) =
      &      hsnow(i,j,bi,bj)*maskC(i,j,1,bi,bj)
+         elseif (gencost_barfile(k)(1:11).EQ.'m_freeboard') then
+           gencost_modfld(i,j,bi,bj,k) =
+     &      gencost_storefld(i,j,bi,bj,k)*maskC(i,j,1,bi,bj)
 #endif
 #ifdef ALLOW_GENCOST3D
          elseif (gencost_barfile(k)(1:7).EQ.'m_theta') then
diff --git a/pkg/ecco/ecco_phys.F b/pkg/ecco/ecco_phys.F
index 94ba6b97cd..b819f55b3a 100644
--- a/pkg/ecco/ecco_phys.F
+++ b/pkg/ecco/ecco_phys.F
@@ -1,4 +1,7 @@
 #include "ECCO_OPTIONS.h"
+#ifdef ALLOW_SEAICE
+# include "SEAICE_OPTIONS.h"
+#endif
 #ifdef ALLOW_SHELFICE
 # include "SHELFICE_OPTIONS.h"
 #endif
@@ -30,6 +33,11 @@ SUBROUTINE ECCO_PHYS( myTime, myIter, myThid )
 # include "ECCO_SIZE.h"
 # include "ECCO.h"
 #endif
+#ifdef ALLOW_SEAICE
+# include "SEAICE_SIZE.h"
+# include "SEAICE.h"
+# include "SEAICE_PARAMS.h"
+#endif
 #ifdef ALLOW_PTRACERS
 # include "PTRACERS_SIZE.h"
 # include "PTRACERS_FIELDS.h"
@@ -82,16 +90,23 @@ SUBROUTINE ECCO_PHYS( myTime, myIter, myThid )
 #ifdef ALLOW_PSBAR_STERIC
       _RL VOLsumTile(nSx,nSy), RHOsumTile(nSx,nSy)
       _RL VOLsumGlob_1, RHOsumGlob_1
-c     CHARACTER*(MAX_LEN_MBUF) msgBuf
+C     CHARACTER*(MAX_LEN_MBUF) msgBuf
 #endif
 C     Mload     :: total mass load (kg/m**2)
-c     _RL Mload(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
+C     _RL Mload(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
+#ifdef ALLOW_SEAICE
+      _RL tmprecip_area,  tmphsnowact, tmpheffact
+      _RL area_reg_sq
+#endif
 #if ( defined ALLOW_AUTODIFF_TAMC && defined ALLOW_GENCOST_CONTRIBUTION )
 C     ikey   :: tape key (depends on kgen)
       INTEGER ikey
 #endif
 CEOP
 
+#ifdef ALLOW_SEAICE
+      area_reg_sq = SEAICE_area_reg * SEAICE_area_reg
+#endif
       tmpfac = recip_rhoConst*recip_gravity
       sIceLoadFacLoc = zeroRL
       IF ( useRealFreshWaterFlux ) sIceLoadFacLoc = recip_rhoConst
@@ -345,6 +360,28 @@ SUBROUTINE ECCO_PHYS( myTime, myIter, myThid )
        ENDDO
        areavolGlob=0. _d 0
 
+#ifdef ALLOW_SEAICE
+       IF (gencost_barfile(kgen)(1:11).EQ.'m_freeboard') THEN
+        DO bj=myByLo(myThid),myByHi(myThid)
+         DO bi=myBxLo(myThid),myBxHi(myThid)
+          DO j = 1,sNy
+           DO i =  1,sNx
+            tmprecip_area = 1. _d 0/SQRT(area(i,j,bi,bj)*area(i,j,bi,bj)
+     &       + area_reg_sq)
+            tmphsnowact = hsnow(I,J,bi,bj) * tmprecip_area
+            tmpheffact  = heff(I,J,bi,bj) * tmprecip_area
+            gencost_storefld(i,j,bi,bj,kgen) =
+     &       ( ( tmphsnowact + tmpheffact ) -
+     &         ( tmphsnowact * SEAICE_rhoSnow +
+     &           tmpheffact * SEAICE_rhoIce ) * recip_rhoConst
+     &       ) * maskC(i,j,1,bi,bj)
+           ENDDO
+          ENDDO
+         ENDDO
+        ENDDO
+       ENDIF
+#endif
+
        DO bj=myByLo(myThid),myByHi(myThid)
         DO bi=myBxLo(myThid),myBxHi(myThid)
          DO j = 1,sNy