From 8ac6ac90df9d5e5eb436b6af5a4f2b333d884551 Mon Sep 17 00:00:00 2001
From: domingom <domingom@ucar.edu>
Date: Wed, 11 Mar 2026 14:57:20 -0600
Subject: [PATCH] Use TKE_0 for boundary conditions in LAD runs

---
 SRC/HYDRO_CORE/CUDA/cuda_hydroCoreDevice.cu | 42 ++++++++++++++-------
 SRC/HYDRO_CORE/hydro_core.c                 | 16 ++++----
 2 files changed, 36 insertions(+), 22 deletions(-)

diff --git a/SRC/HYDRO_CORE/CUDA/cuda_hydroCoreDevice.cu b/SRC/HYDRO_CORE/CUDA/cuda_hydroCoreDevice.cu
index 7c85ccd..3fc8934 100644
--- a/SRC/HYDRO_CORE/CUDA/cuda_hydroCoreDevice.cu
+++ b/SRC/HYDRO_CORE/CUDA/cuda_hydroCoreDevice.cu
@@ -761,18 +761,34 @@ __global__ void cudaDevice_hydroCoreCommence(int simTime_it, float* hydroFlds_d,
        fldBS = &sgstkeScalarsFrhs_d[fldStride*iFld]; // Frhs forcing iwas set to zero, so it can be used here as zero-valued base state
        if(hydroBCs_d == 1){ //Using LAD BCs
 	cudaDevice_VerticalAblZeroGradBCs(fld);
-	if(rankXid_d == 0){
+	if (iFld == 0){ // TKE_0
+	  if(rankXid_d == 0){
+           cudaDevice_westBdyBCs(iFld+Nhydro_d, timeWeight, fld, YZBdyPlanes_d, YZBdyPlanesNext_d);
+          }
+          if(rankXid_d == numProcsX_d-1){
+           cudaDevice_eastBdyBCs(iFld+Nhydro_d, timeWeight, fld, YZBdyPlanes_d, YZBdyPlanesNext_d);
+          }
+          if(rankYid_d == 0){
+           cudaDevice_southBdyBCs(iFld+Nhydro_d, timeWeight, fld, XZBdyPlanes_d, XZBdyPlanesNext_d);
+          }
+          if(rankYid_d == numProcsY_d-1){
+           cudaDevice_northBdyBCs(iFld+Nhydro_d, timeWeight, fld, XZBdyPlanes_d, XZBdyPlanesNext_d);
+          }
+          cudaDevice_ceilingBdyBCs(iFld+Nhydro_d, timeWeight, fld, XYBdyPlanes_d, XYBdyPlanesNext_d);
+	}else{ // all other TKE scales
+	  if(rankXid_d == 0){
            cudaDevice_lateralTKEBdyBCs(iFld, fld, fldBS, 0);
-         }
-         if(rankXid_d == numProcsX_d-1){
+          }
+          if(rankXid_d == numProcsX_d-1){
            cudaDevice_lateralTKEBdyBCs(iFld, fld, fldBS, 1);
-         }
-         if(rankYid_d == 0){
+          }
+          if(rankYid_d == 0){
            cudaDevice_lateralTKEBdyBCs(iFld, fld, fldBS, 2);
-         }
-         if(rankYid_d == numProcsY_d-1){
+          }
+          if(rankYid_d == numProcsY_d-1){
            cudaDevice_lateralTKEBdyBCs(iFld, fld, fldBS, 3);
-         }
+          }
+	} // end if (iFld == 0)
        }else if (hydroBCs_d == 2){
 	 cudaDevice_VerticalAblZeroGradBCs(fld);
          if(numProcsX_d==1){
@@ -795,18 +811,18 @@ __global__ void cudaDevice_hydroCoreCommence(int simTime_it, float* hydroFlds_d,
        if(hydroBCs_d == 1){ //Using LAD BCs
          cudaDevice_VerticalAblBCs(iFld, fld, fldBS);
          if(rankXid_d == 0){
-           cudaDevice_westBdyBCs(iFld+Nhydro_d, timeWeight, fld, YZBdyPlanes_d, YZBdyPlanesNext_d);
+           cudaDevice_westBdyBCs(iFld+Nhydro_d+1, timeWeight, fld, YZBdyPlanes_d, YZBdyPlanesNext_d);
          }
          if(rankXid_d == numProcsX_d-1){
-           cudaDevice_eastBdyBCs(iFld+Nhydro_d, timeWeight, fld, YZBdyPlanes_d, YZBdyPlanesNext_d);
+           cudaDevice_eastBdyBCs(iFld+Nhydro_d+1, timeWeight, fld, YZBdyPlanes_d, YZBdyPlanesNext_d);
          }
          if(rankYid_d == 0){
-           cudaDevice_southBdyBCs(iFld+Nhydro_d, timeWeight, fld, XZBdyPlanes_d, XZBdyPlanesNext_d);
+           cudaDevice_southBdyBCs(iFld+Nhydro_d+1, timeWeight, fld, XZBdyPlanes_d, XZBdyPlanesNext_d);
          }
          if(rankYid_d == numProcsY_d-1){
-           cudaDevice_northBdyBCs(iFld+Nhydro_d, timeWeight, fld, XZBdyPlanes_d, XZBdyPlanesNext_d);
+           cudaDevice_northBdyBCs(iFld+Nhydro_d+1, timeWeight, fld, XZBdyPlanes_d, XZBdyPlanesNext_d);
          }
-         cudaDevice_ceilingBdyBCs(iFld+Nhydro_d, timeWeight, fld, XYBdyPlanes_d, XYBdyPlanesNext_d);
+         cudaDevice_ceilingBdyBCs(iFld+Nhydro_d+1, timeWeight, fld, XYBdyPlanes_d, XYBdyPlanesNext_d);
        }else if (hydroBCs_d == 2){
          cudaDevice_VerticalAblZeroGradBCs(fld); // to apply zero-gradient bottom/top BCs
          if(numProcsX_d==1){
diff --git a/SRC/HYDRO_CORE/hydro_core.c b/SRC/HYDRO_CORE/hydro_core.c
index 2bbefe6..25d6e3a 100644
--- a/SRC/HYDRO_CORE/hydro_core.c
+++ b/SRC/HYDRO_CORE/hydro_core.c
@@ -1551,12 +1551,12 @@ int hydro_coreInit(){
        fflush(stdout);
      }
      if( moistureSelector > 0){
-       nBndyVars = Nhydro+moistureNvars;
+       nBndyVars = Nhydro+1+moistureNvars; // +1 is for TKE_0
        nSurfBndyVars = 2;   //Only allows tskin and qskin
      }else{
-       nBndyVars = Nhydro;
+       nBndyVars = Nhydro+1; // +1 is for TKE_0
        nSurfBndyVars = 1;   //Only allows tskin
-     } //end if moisture is on else not   //NOTE: Doesn't handle any AuxScalars or TKE-related Prog. variables.
+     } //end if moisture is on else not   //NOTE: Doesn't handle any AuxScalars Prog. variables.
      XZBdyPlanesGlobal = (float *) malloc( 2*(nBndyVars)*Nx*Nz*sizeof(float) );
      YZBdyPlanesGlobal = (float *) malloc( 2*(nBndyVars)*Ny*Nz*sizeof(float) );
      XYBdyPlanesGlobal = (float *) malloc( 2*(nBndyVars)*Nx*Ny*sizeof(float) );
@@ -2000,19 +2000,17 @@ int hydro_coreSetupBndyPlanesAllRanks(){
        sprintf(fieldName,"theta");
        fieldIndex = 4;
        errorCode = hydro_coreReadFieldBndyPlanes(ncid, fieldName, fieldIndex);
+       sprintf(fieldName,"TKE_0");
+       fieldIndex = 5;
+       errorCode = hydro_coreReadFieldBndyPlanes(ncid, fieldName, fieldIndex);
        if(moistureSelector > 0){
          if(moistureNvars > 0){
            sprintf(fieldName,"qv");
-           fieldIndex = 5;
+           fieldIndex = 6;
            errorCode = hydro_coreReadFieldBndyPlanes(ncid, fieldName, fieldIndex);
          }
          if(moistureNvars > 1){
            sprintf(fieldName,"ql");
-           fieldIndex = 6;
-           errorCode = hydro_coreReadFieldBndyPlanes(ncid, fieldName, fieldIndex);
-         }
-         if(moistureNvars > 2){
-           sprintf(fieldName,"qr");
            fieldIndex = 7;
            errorCode = hydro_coreReadFieldBndyPlanes(ncid, fieldName, fieldIndex);
          }