fixup compute_land_ice_draft_from_pressure

Author: cbegeman

compass/ocean/iceshelf.py

@@ -96,10 +96,8 @@ def compute_land_ice_draft_from_pressure(land_ice_pressure, modify_mask,
     if ref_density is None:
         ref_density = constants['SHR_CONST_RHOSW']
     land_ice_draft_array = \
-        - (modify_mask.values *
-           land_ice_pressure.values / (ref_density * gravity))
-    land_ice_draft = xarray.DataArray(data=land_ice_draft_array,
-                                      dims=(land_ice_pressure.dims))
+        modify_mask * -land_ice_pressure / (ref_density * gravity)
+    land_ice_draft = land_ice_draft_array.transpose(*land_ice_pressure.dims)
     return land_ice_draft
 
 

compass/ocean/tests/isomip_plus/initial_state.py

@@ -158,7 +158,8 @@ def _compute_initial_condition(self):
             modify_mask = ds.bottomDepth > 0.
             land_ice_draft = compute_land_ice_draft_from_pressure(
                 land_ice_pressure=land_ice_pressure,
-                modify_mask=modify_mask)
+                modify_mask=modify_mask,
+                ref_density=1028.)
             land_ice_draft = np.maximum(land_ice_draft,
                                         -ds.bottomDepth)
             ds['ssh'] = land_ice_draft