@@ -167,45 +167,49 @@ def get_model(ws, name, array_input=False):
167167 drnspd = {}
168168 if array_input :
169169 bheadspd = {}
170- condspd = {}
171- auxspd = {}
170+ ghbcondspd = {}
171+ ghbauxspd = {}
172172 else :
173173 ghbspd = {}
174174 for kper in range (nper ):
175175 if kper == 0 :
176176 sl = sealevel
177177 else :
178178 sl = sealevelts [kper ]
179+ sl = np .round (sl , decimals = 8 )
179180 drnlist = []
180181 if array_input :
181- abhead = np .full ((nlay , nrow , ncol ), DNODATA , dtype = float )
182- acond = np .full ((nlay , nrow , ncol ), DNODATA , dtype = float )
183- aconc = np .full ((nlay , nrow , ncol ), DNODATA , dtype = float )
184- adens = np .full ((nlay , nrow , ncol ), DNODATA , dtype = float )
182+ bhead = np .full ((nlay , nrow , ncol ), DNODATA , dtype = float )
183+ ghbcond = np .full ((nlay , nrow , ncol ), DNODATA , dtype = float )
184+ ghbconc = np .full ((nlay , nrow , ncol ), DNODATA , dtype = float )
185+ ghbdens = np .full ((nlay , nrow , ncol ), DNODATA , dtype = float )
185186 else :
186187 ghblist = []
187- nbound = 0
188+ ghbbnd = 0
189+ drnbnd = 0
188190 for k , i , j in zip (kidx , iidx , jidx ):
189191 zcell = zcellcenters [k , i , j ]
190192 cond = 864.0 * (delz * delc ) / (0.5 * delr )
191193 if zcell > sl :
192194 drnlist .append ([(k , i , j ), zcell , 864.0 , 0.0 ])
195+ drnbnd += 1
193196 else :
194197 if array_input :
195- abhead [k , i , j ] = sl
196- acond [k , i , j ] = 864.0
197- aconc [k , i , j ] = 35.0
198- adens [k , i , j ] = 1024.5
198+ bhead [k , i , j ] = sl
199+ ghbcond [k , i , j ] = 864.0
200+ ghbconc [k , i , j ] = 35.0
201+ ghbdens [k , i , j ] = 1024.5
199202 else :
200203 ghblist .append ([(k , i , j ), sl , 864.0 , 35.0 , 1024.5 ])
201- nbound += 1
202- if array_input and zcell <= sl :
203- bheadspd [kper ] = abhead
204- condspd [kper ] = acond
205- auxspd [kper ] = [aconc , adens ]
206- elif len (ghblist ) > 0 :
207- ghbspd [kper ] = ghblist
208- if len (drnlist ) > 0 :
204+ ghbbnd += 1
205+ if ghbbnd > 0 :
206+ if array_input :
207+ bheadspd [kper ] = bhead
208+ ghbcondspd [kper ] = ghbcond
209+ ghbauxspd [kper ] = [ghbconc , ghbdens ]
210+ else :
211+ ghbspd [kper ] = ghblist
212+ if drnbnd > 0 :
209213 drnspd [kper ] = drnlist
210214
211215 # drn
@@ -230,8 +234,8 @@ def get_model(ws, name, array_input=False):
230234 pname = "GHB-1" ,
231235 auxiliary = ["CONCENTRATION" , "DENSITY" ],
232236 bhead = bheadspd ,
233- cond = condspd ,
234- aux = auxspd ,
237+ cond = ghbcondspd ,
238+ aux = ghbauxspd ,
235239 )
236240 else :
237241 ghb1 = flopy .mf6 .ModflowGwfghb (
0 commit comments