From aed27ca548059dfe41a175a23b3a93cb654be8fe Mon Sep 17 00:00:00 2001
From: Alice Barthel <abarthel@lanl.gov>
Date: Mon, 6 Nov 2023 19:27:07 -0600
Subject: [PATCH] added config parameters for initial temp profile

---
 .../tests/baroclinic_gyre/baroclinic_gyre.cfg | 12 +++++---
 .../tests/baroclinic_gyre/initial_state.py    | 30 ++++++++++++++-----
 2 files changed, 30 insertions(+), 12 deletions(-)

diff --git a/compass/ocean/tests/baroclinic_gyre/baroclinic_gyre.cfg b/compass/ocean/tests/baroclinic_gyre/baroclinic_gyre.cfg
index d3187cad23..54c2364b9c 100644
--- a/compass/ocean/tests/baroclinic_gyre/baroclinic_gyre.cfg
+++ b/compass/ocean/tests/baroclinic_gyre/baroclinic_gyre.cfg
@@ -32,15 +32,19 @@ lat_max = 75
 lon_min = 0
 lon_max = 60
 
-# Maximum zonal wind stress value
-wind_stress_max = 0.1
+# Initial vertical temperature profile
+initial_temp_top = 33.
+initial_temp_bot = 1.
 
 # Constant salinity value (also used in restoring)
 initial_salinity = 34.
 
+# Maximum zonal wind stress value
+wind_stress_max = 0.1
+
 # Surface temperature restoring
-temp_min = 0.
-temp_max = 30.
+restoring_temp_min = 0.
+restoring_temp_max = 30.
 
 # Restoring piston velocity for surface temperature (s-1)
 restoring_temp_piston_vel = 3.86e-7
diff --git a/compass/ocean/tests/baroclinic_gyre/initial_state.py b/compass/ocean/tests/baroclinic_gyre/initial_state.py
index 9f72ae74f1..23b2f85382 100644
--- a/compass/ocean/tests/baroclinic_gyre/initial_state.py
+++ b/compass/ocean/tests/baroclinic_gyre/initial_state.py
@@ -54,7 +54,6 @@ def run(self):
 
 
 def _write_initial_state(config, dsMesh):
-    # section = config['baroclinic_gyre']
 
     ds = dsMesh.copy()
 
@@ -65,12 +64,27 @@ def _write_initial_state(config, dsMesh):
     init_vertical_coord(config, ds)
 
     # setting the initial conditions
-    initial_salinity = config.getfloat('baroclinic_gyre', 'initial_salinity')
-    temperature = (-11. * np.log(0.0414 *
-                   (-1. * ds.zMid + 100.3)) + 48.8)
+    section = config['baroclinic_gyre']
+    initial_salinity = section.getfloat('initial_salinity')
+    temp_top = section.getfloat('initial_temp_top')
+    temp_bot = section.getfloat('initial_temp_bot')
+
+    cc = 0.049  # 0.049 value from optimization on 1800m
+    aa = temp_top / np.log(cc)
+    bb = (cc ** (temp_bot / temp_top) - cc)
+    temperature = aa * np.log(bb * ds.zMid / bottom_depth + cc)
+    val_bot = aa * np.log(bb + cc)
+    val_top = aa * np.log(cc)
+    print(f'analytical bottom T: {val_bot} at \
+            depth : {bottom_depth}')
+    print(f'analytical surface T: {val_top} at \
+            depth = 0')
+
+#    temperature = (-11. * np.log(0.0414 *
+#                   (-1. * ds.zMid + 100.3)) + 48.8)
     temperature = temperature.transpose('Time', 'nCells', 'nVertLevels')
-    print(f'bottom T: {temperature[0, 0, -1]} and \
-            surface : {temperature[0, 0, 0]}')
+    print(f'bottom layer T: {temperature[0, 0, -1]} and \
+            surface layer T: {temperature[0, 0, 0]}')
     salinity = initial_salinity * xr.ones_like(temperature)
 
     normalVelocity = xr.zeros_like(ds.xEdge)
@@ -96,8 +110,8 @@ def _write_forcing(config, lat):
     latMin = section.getfloat('lat_min')
     latMax = section.getfloat('lat_max')
     tauMax = section.getfloat('wind_stress_max')
-    tempMin = section.getfloat('temp_min')
-    tempMax = section.getfloat('temp_max')
+    tempMin = section.getfloat('restoring_temp_min')
+    tempMax = section.getfloat('restoring_temp_max')
     restoring_temp_piston_vel = section.getfloat('restoring_temp_piston_vel')
     initial_salinity = section.getfloat('initial_salinity')
     lat = np.rad2deg(lat)