From b0b72654ea5b04a900def9e43318b8ad1ac7cad6 Mon Sep 17 00:00:00 2001
From: Bruce Perry <bruce.perry@nrel.gov>
Date: Tue, 16 Jan 2024 18:13:56 -0700
Subject: [PATCH] terrain initialization for GABLS1

---
 Exec/RegTests/GABLS1/prob.cpp | 19 ++++++++++++++++---
 1 file changed, 16 insertions(+), 3 deletions(-)

diff --git a/Exec/RegTests/GABLS1/prob.cpp b/Exec/RegTests/GABLS1/prob.cpp
index 5aed34a33..436da69ce 100644
--- a/Exec/RegTests/GABLS1/prob.cpp
+++ b/Exec/RegTests/GABLS1/prob.cpp
@@ -55,7 +55,7 @@ Problem::init_custom_pert(
     amrex::Array4<amrex::Real      > const& r_hse,
     amrex::Array4<amrex::Real      > const& /*p_hse*/,
     amrex::Array4<amrex::Real const> const& /*z_nd*/,
-    amrex::Array4<amrex::Real const> const& /*z_cc*/,
+    amrex::Array4<amrex::Real const> const& z_cc,
     amrex::GeometryData const& geomdata,
     amrex::Array4<amrex::Real const> const& /*mf_m*/,
     amrex::Array4<amrex::Real const> const& /*mf_u*/,
@@ -71,21 +71,35 @@ Problem::init_custom_pert(
     amrex::Print() << "Initializing QKE" << std::endl;
   }
 
+  const bool use_terrain = sc.use_terrain;
+
   if (parms.pert_ref_height > 0) {
     if (parms.pert_deltaU != 0.0) {
       amrex::Print() << "Adding divergence-free x-velocity perturbations" << std::endl;
+      if (use_terrain) {
+          amrex::Abort("u perturbations not supported for terrain");
+      }
     }
     if (parms.pert_deltaV != 0.0) {
       amrex::Print() << "Adding divergence-free y-velocity perturbations" << std::endl;
+      if (use_terrain) {
+          amrex::Abort("v perturbations not supported for terrain");
+      }
     }
     if (parms.T_0_Pert_Mag != 0.0) {
       amrex::Print() << "Adding random theta perturbations" << std::endl;
     }
     if (parms.U_0_Pert_Mag != 0.0) {
       amrex::Print() << "Adding random x-velocity perturbations" << std::endl;
+      if (use_terrain) {
+          amrex::Abort("u perturbations not supported for terrain");
+      }
     }
     if (parms.V_0_Pert_Mag != 0.0) {
       amrex::Print() << "Adding random y-velocity perturbations" << std::endl;
+      if (use_terrain) {
+          amrex::Abort("v perturbations not supported for terrain");
+      }
     }
   }
 
@@ -96,7 +110,7 @@ Problem::init_custom_pert(
     const Real* dx = geomdata.CellSize();
     const Real x = prob_lo[0] + (i + 0.5) * dx[0];
     const Real y = prob_lo[1] + (j + 0.5) * dx[1];
-    const Real z = prob_lo[2] + (k + 0.5) * dx[2];
+    const Real z = use_terrain ? z_cc(i,j,k) : prob_lo[2] + (k + 0.5) * dx[2];
 
     // Define a point (xc,yc,zc) at the center of the domain
     const Real xc = 0.5 * (prob_lo[0] + prob_hi[0]);
@@ -196,4 +210,3 @@ Problem::init_custom_pert(
     }
   });
 }
-