added config parameters for initial temp profile

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

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)