# Python script linking SWAT, MODFLOW and MT3D models
# for the SOILPROM project Use case 8: Puck Bay
# version under development

import csv
import numpy as np
import flopy
import matplotlib.pyplot as plt
import flopy.utils as fu
import flopy.utils.binaryfile as bf
import flopy.utils.formattedfile as ff
import subprocess

# names of the executable files if each model

swat_exe_name = 'SWAT2012.exe'
modflow_exe_name = 'MODFLOW-NWT_64.exe'
mt3dms_exe_name = 'mt3dms5b.exe'

subprocess.run([swat_exe_name])

# spatial coordinates of the model domain for MODFLOW (rectangle)
xNW = 451730
yNW = 770900

xSE = 473780
ySE = 747950

# model grid has 459 rows and 441 columns
# square gridblocks 50x50 m

# reading MODFLOW model information
ml = flopy.modflow.Modflow.load('WP.nam',verbose=True)

my_nlay = ml.dis.nlay
my_nrow = ml.dis.nrow
my_ncol = ml.dis.ncol
my_nper = ml.dis.nper

# reading initial condition for MODFLOW
# (hydraulic head) from file WP.hds
# after the simulation the file will be overwritten
# with values calculated for the current day
# and can be used as the initial condition for the next day 
hfile = bf.HeadFile('WP.hds')
# get values from the last time step
end_time = hfile.get_times()[-1]
init_head = hfile.get_data(totim=end_time)
# copy existing array IBOUND
my_ibound = ml.bas6.ibound.array.copy()
# create new array STRT
my_strt = init_head

# modify BAS package - new initial condition
bas6 = flopy.modflow.mfbas.ModflowBas(model=ml,ibound=my_ibound,strt=my_strt)
ml.bas6.write_file()


# read recharge date from file output.hru generated by SWAT
# for a specific day, e.g. 1 Feb 2001
my_day = 1
my_month = 2
my_year = 2001
# number of HRUs in SWAT model (constant)
n_hru = 564
# initialize array with values for HRUs
swat_data = np.zeros((n_hru,6))
with open ('output.hru') as out_file:
  out_reader = csv.reader(out_file, delimiter=' ', skipinitialspace=True)
  for i in range(9):  # skip the heading
    next(out_reader,None)
  while True:
    row = next(out_reader,None)
    if row == None:
      break
    r_hru = int(row[1])-1
    r_month = int(row[5])
    r_day = int(row[6])
    r_year = int(row[7])
    #print(r_day, r_month, r_year)
    r_area = float(row[8]) # area of HRU [km2]
    r_perc = float(row[17]) # recharge [mm/d]
    r_evap = float(row[20]) # ET from groundwater [mm/d]
    r_no3 = float(row[64]) # load N-NO3 [kg/ha/d]
    # check the date
    if (r_day == my_day) and (r_month == my_month) and (r_year == my_year):
      swat_data[r_hru,0]=r_area
      swat_data[r_hru,1]=r_perc
      swat_data[r_hru,2]=r_evap
      swat_data[r_hru,3]=r_no3
      swat_data[r_hru,4]=0
      swat_data[r_hru,5]=0
      if r_perc > r_evap:
        # net recharge [m/h]
        swat_data[r_hru,4]=(r_perc-r_evap)/24000
        # N-NO3 concentration in recharge [mg/dm3]
        swat_data[r_hru,5]=r_no3/swat_data[r_hru,4]/240



# calculate average values for SWAT model area
area_total = np.sum(swat_data[:,0])
rch_aver = np.sum(np.multiply(swat_data[:,0],swat_data[:,4]))/area_total
no3_aver = np.sum(np.multiply(swat_data[:,0],swat_data[:,3]))/area_total
no3_conc_aver = 0
if rch_aver > 0:
  no3_conc_aver = no3_aver/rch_aver/240
print('SWAT model area: ',str(area_total),' km2')
print('Average recharge: ',str(rch_aver), 'm/h')
print('Average N-NO3 conc. in recharge: ',str(no3_conc_aver), 'mg/dm3')

# read files mapping SWAT HRUs to MODFLOW grid cells
# hru_dhru.txt and grid_dhru.txt

# reading file hru_dhru.txt
# indices of HRUs corresponding to specific DHRUs

idhh = []
with open ('hru_dhru.txt') as out_file:
  out_reader = csv.reader(out_file, delimiter='\t', skipinitialspace=True)
  # read number of DHRUs (divides HRU)
  row = next(out_reader,None)
  n_dhru = int(row[0])
  # read number of HRUs
  row = next(out_reader,None)
  #my_hru = int(row[0])
  # read heading
  row = next(out_reader,None)
  for i in range(n_dhru):
    row = next(out_reader,None)
    idhh.append(int(row[2]))

grid_data = []

with open ('grid_dhru.txt') as out_file:
  out_reader = csv.reader(out_file, delimiter='\t', skipinitialspace=True)
  # read the number of records
  row = next(out_reader,None)
  n_grid = int(row[0])
  # read the number of DHRUs (divided HRUs)
  row = next(out_reader,None)
  #n_dhru = int(row[0])
  # read the number of rows
  row = next(out_reader,None) 
  #my_row = int(row[0])
  # read the number of columns
  row = next(out_reader,None)
  #my_ncol = int(row[0])
  # read heading
  row = next(out_reader,None)
  for i in range(n_grid):
    row = next(out_reader,None)
    grid_data.append(row)

#my_grid = int(grid_data[0][0])
#print(my_grid)
#print(ml.dis.get_lrc(my_grid)[0][1])
# preparing array with recharge values
# in all grid cells

my_mask = ml.bas6.ibound.array[2,:,:].copy()
my_rch = np.zeros((my_nrow,my_ncol))
# recharge equal to the average value for the whole area
my_rch[my_mask>0] = rch_aver
my_no3_load = np.zeros((my_nrow,my_ncol))
# nitrate load equal to the average value for the whole area
my_no3_load[my_mask>0] = no3_conc_aver*rch_aver

for igrid in grid_data:
  my_grid = int(igrid[0])
# find the row and column index of the given grid cell
  my_row = ml.dis.get_lrc(my_grid)[0][1]
  my_col = ml.dis.get_lrc(my_grid)[0][2]
  #print(my_grid,my_row,my_col)
# check if the cell is active
  if my_mask[my_row,my_col] > 0:
    # modify recharge and concentration
    grid_area = float(igrid[1])
    overlap = float(igrid[3])
    coeff = overlap/grid_area
    #print(grid_area,overlap,coeff)
    # find DHRU index, indexing from 0
    my_dhru = int(igrid[2])-1
    # find HRU index, indexing from 0
    my_hru = idhh[my_dhru]-1
    #print(my_dhru,my_hru)
    hru_rch = swat_data[my_hru,4]
    hru_cno3 = swat_data[my_hru,5]
    new_rch = my_rch[my_row,my_col]+(hru_rch-rch_aver)*coeff
    new_no3_load = (my_no3_load[my_row,my_col]\
               +(hru_cno3*hru_rch-rch_aver*no3_conc_aver)*coeff)
    #print(my_rch[my_row,my_col])
    #print(hru_rch,new_rch)
    my_rch[my_row,my_col] = new_rch
    my_no3_load[my_row,my_col] = new_no3_load
    #print(my_rch[my_row,my_col])
    #input('Press Enter')
my_no3_conc = np.zeros((my_nrow,my_ncol))
my_no3_conc[my_rch>0] = np.divide(my_no3_load[my_rch>0],my_rch[my_rch>0])


# adding new recharge values to the model
flopy.modflow.mfrch.ModflowRch(ml, nrchop=3, ipakcb=40, \
    rech=my_rch, unitnumber=10)
ml.rch.write_file()


ml.exe_name = modflow_exe_name
# running MODFLOW program
success, buff = ml.run_model()
print(success)
print(buff)



# reading MT3D model
mt = flopy.mt3d.mt.Mt3dms.load('WP.mt_nam',exe_name=mt3dms_exe_name, verbose=True)
# specify MODFLOW file which is the basis of the transport simulation
mt.ftlfilename = 'WP.ftl'

# read N-NO3 concentration (from the previous step)
# file WP_NO3.ucn, which will be later overwritten with the results for the current time step
my_c0 = np.zeros((my_nlay, my_nrow, my_ncol))
cfile = bf.UcnFile('WP_NO3.ucn')
# control listing of time steps
#for t in cfile.get_times():
#  print(str(t))
# results from the last time step
end_time = cfile.get_times()[-1]
#print(str(end_time))
for ilay in range(my_nlay):
  my_c0[ilay,:,:] = cfile.get_data(totim=end_time,mflay=ilay)
#cN_2 = cfile.get_data(totim=end_time,mflay=2)
#cN_4 = cfile.get_data(totim=end_time,mflay=4)

# writing new initial condition for MT3D (concentrations) to file WP.btn

mt.btn.sconc = []
u3d = flopy.utils.util_array.Util3d(mt, (my_nlay, my_nrow, my_ncol), np.float32,
                     my_c0, name='sconc1', locat=mt.btn.unit_number[0],
                     array_free_format=False)
mt.btn.sconc.append(u3d)
mt.btn.write_file()

# writing new values of N-NO3 load to file WP.ssm

mt.ssm.crch = []
t2d = flopy.utils.util_array.Transient2d(mt, (my_nrow, my_ncol), np.float32, 
                    my_no3_conc, name='crch1',
                    locat=mt.ssm.unit_number[0], array_free_format=False)
mt.ssm.crch.append(t2d)
mt.ssm.write_file()

# modify the first row in WP.ssm, 
# necessary to read the file in subsequent steps
ssm_file = open('WP.ssm','rt') # open file
ssm_lines = ssm_file.readlines() # read all lines
ssm_file.close() # close the file
ssm_lines[0]=' F F T F F T F F F F\n'
ssm_file = open('WP.ssm','wt') # open file for writing
for line in ssm_lines: # write line
  ssm_file.write(line)
ssm_file.close() # close the file

# running MT3D computer program
success, buff = mt.run_model()
print(success)
print(buff)

# reading fluxes between grid cells
ffile = bf.CellBudgetFile('WP.cbc')
#ffile.list_records()
    
# control listing of time steps
#for t in ffile.get_times():
#  print(str(t))
# results from the last time step
end_time = ffile.get_times()[0]
#print(str(end_time))
# flow in x firection (W-E)
flow1 = ffile.get_data(text='FLOW RIGHT FACE', totim=end_time,full3D=True)
# flow in y direction (N-S)
flow2 = ffile.get_data(text='FLOW FRONT FACE', totim=end_time,full3D=True)
# flow in z direction
flow3 = ffile.get_data(text='FLOW LOWER FACE', totim=end_time,full3D=True)

# read the calculated concentration of N-NO3
# we are interested in layers "2" and "4" (indexing from 0)
# kontrolny wydruk kroków czasowych z pliku wynikowego
cfile = bf.UcnFile('WP_NO3.ucn')
# control listing of time steps
#for t in cfile.get_times():
#  print(str(t))
# results from the last time step 
end_time = cfile.get_times()[0]
#print(str(end_time))
cN_2 = cfile.get_data(totim=end_time,mflay=2)
cN_4 = cfile.get_data(totim=end_time,mflay=4)


# control plots of hydraulic head and concentration

hdata = bf.HeadFile('WP.hds')
# get data from the last time step
end_time = hdata.get_times()[-1]
my_heads = hdata.get_data(totim=end_time)

#fig = plt.figure(figsize=(10, 6))
#ax = fig.add_subplot(1, 2, 1, aspect="equal")
#modelmap = flopy.plot.PlotMapView(model=ml, layer=0, ax=ax)
#pa1 = modelmap.plot_array(new_head, masked_values=[-999.], cmap='rainbow')
#cb = plt.colorbar(pa1, shrink=0.5)

plot_heads = np.copy(my_heads[2, :, :])
plot_heads[my_heads[2, :, :] < -900.0] = np.nan
plot1 = plt.figure(1)
plt.subplot(1, 1, 1, aspect='equal')
plt.title('Head distribution (m)')
plt.imshow(plot_heads, cmap='rainbow', vmin = 0.0)
plt.colorbar()

#contours = plt.contour(np.flipud(head[0, :, :]), levels=levels, extent=extent, zorder=10)
#plt.clabel(contours, inline=1, fontsize=10, fmt='%d', zorder=11)

#plt.show()

plot2 = plt.figure(2)
cdata = fu.UcnFile('WP_NO3.UCN')

end_time = cdata.get_times()[-1]
conc = cdata.get_data(totim=end_time, mflay=2)
conc[my_heads[2, :, :] < -900.0] = np.nan
#conc.plot(totim=end_time, colorbar='NO3 concentration (mg/l)', cmap='Blues', vmin=0.)
plt.imshow(conc, cmap='rainbow', vmin = 0.0)
plt.title('Concentration distribution (mg/l)')
plt.colorbar()
plt.show()


