# this runs CLM test case

#
# Import the ParFlow TCL package
#
lappend auto_path $env(PARFLOW_DIR)/bin 
package require parflow
namespace import Parflow::*

foreach dir {qflx_evap_grnd eflx_lh_tot qflx_evap_tot qflx_tran_veg correct_output qflx_infl swe_out eflx_lwrad_out t_grnd diag_out qflx_evap_soi eflx_soil_grnd eflx_sh_tot qflx_evap_veg qflx_top_soil} {
    file mkdir $dir
}

#-----------------------------------------------------------------------------
# File input version number
#-----------------------------------------------------------------------------
pfset FileVersion 4

#-----------------------------------------------------------------------------
# Process Topology
#-----------------------------------------------------------------------------

set P  [lindex $argv 0]
set Q  [lindex $argv 1]
set R  [lindex $argv 2]

set NumPatches [lindex $argv 3]

pfset Process.Topology.P $P
pfset Process.Topology.Q $Q   
pfset Process.Topology.R $R

set NumProcs [expr $P * $Q * $R]

#-----------------------------------------------------------------------------
# Computational Grid
#-----------------------------------------------------------------------------

pfset ComputationalGrid.Lower.X                0.0
pfset ComputationalGrid.Lower.Y                0.0
pfset ComputationalGrid.Lower.Z                 0.0

pfset ComputationalGrid.DX	               1000.
pfset ComputationalGrid.DY                     1000. 
pfset ComputationalGrid.DZ	                 0.5

pfset ComputationalGrid.NX                      5
pfset ComputationalGrid.NY                      5
pfset ComputationalGrid.NZ                     10 

if {[expr $NumProcs == 1]} {
    if {[expr $NumPatches == 1]} {
	pfset ProcessGrid.NumSubgrids 1
	pfset ProcessGrid.0.P 0
	pfset ProcessGrid.0.IX 0
	pfset ProcessGrid.0.IY 0
	pfset ProcessGrid.0.IZ 0
	
	pfset ProcessGrid.0.NX 5
	pfset ProcessGrid.0.NY 5
	pfset ProcessGrid.0.NZ 10
    } elseif {[expr $NumPatches == 2]} {
	pfset ProcessGrid.NumSubgrids 2
	pfset ProcessGrid.0.P 0
	pfset ProcessGrid.0.IX 0
	pfset ProcessGrid.0.IY 0
	pfset ProcessGrid.0.IZ 0
	
	pfset ProcessGrid.0.NX 2
	pfset ProcessGrid.0.NY 5
	pfset ProcessGrid.0.NZ 10
	
	pfset ProcessGrid.1.P 0
	pfset ProcessGrid.1.IX 2
	pfset ProcessGrid.1.IY 0
	pfset ProcessGrid.1.IZ 0
	
	pfset ProcessGrid.1.NX 3
	pfset ProcessGrid.1.NY 5
	pfset ProcessGrid.1.NZ 10
    } elseif {[expr $NumPatches == 3]} {
	
	pfset ProcessGrid.NumSubgrids 3
	pfset ProcessGrid.0.P 0
	pfset ProcessGrid.0.IX 0
	pfset ProcessGrid.0.IY 0
	pfset ProcessGrid.0.IZ 0
	pfset ProcessGrid.0.NX 5
	pfset ProcessGrid.0.NY 5
	pfset ProcessGrid.0.NZ 8
	
	pfset ProcessGrid.1.P 0
	pfset ProcessGrid.1.IX 0
	pfset ProcessGrid.1.IY 5
	pfset ProcessGrid.1.IZ 0
	pfset ProcessGrid.1.NX 5
	pfset ProcessGrid.1.NY 5
	pfset ProcessGrid.1.NZ 8

	pfset ProcessGrid.2.P 0
	pfset ProcessGrid.2.IX 5
	pfset ProcessGrid.2.IY 0
	pfset ProcessGrid.2.IZ 0
	pfset ProcessGrid.2.NX 5
	pfset ProcessGrid.2.NY 10
	pfset ProcessGrid.2.NZ 8
    } else {
	puts "Invalid processor/number of subgrid option"
	exit
    }
}


#-----------------------------------------------------------------------------
# The Names of the GeomInputs
#-----------------------------------------------------------------------------
pfset GeomInput.Names "domain_input"


#-----------------------------------------------------------------------------
# Domain Geometry Input
#-----------------------------------------------------------------------------
pfset GeomInput.domain_input.InputType            Box
pfset GeomInput.domain_input.GeomName             domain

#-----------------------------------------------------------------------------
# Domain Geometry
#-----------------------------------------------------------------------------
pfset Geom.domain.Lower.X                        0.0 
pfset Geom.domain.Lower.Y                        0.0
pfset Geom.domain.Lower.Z                          0.0

pfset Geom.domain.Upper.X                        5000.
pfset Geom.domain.Upper.Y                        5000.
pfset Geom.domain.Upper.Z                       5. 

pfset Geom.domain.Patches  "x-lower x-upper y-lower y-upper z-lower z-upper"

#-----------------------------------------------------------------------------
# Perm
#-----------------------------------------------------------------------------
pfset Geom.Perm.Names "domain"

pfset Geom.domain.Perm.Type            Constant
pfset Geom.domain.Perm.Value           0.2


pfset Perm.TensorType               TensorByGeom

pfset Geom.Perm.TensorByGeom.Names  "domain"

pfset Geom.domain.Perm.TensorValX  1.0
pfset Geom.domain.Perm.TensorValY  1.0
pfset Geom.domain.Perm.TensorValZ  1.0

#-----------------------------------------------------------------------------
# Specific Storage
#-----------------------------------------------------------------------------
# specific storage does not figure into the impes (fully sat) case but we still
# need a key for it

pfset SpecificStorage.Type            Constant
pfset SpecificStorage.GeomNames       "domain"
pfset Geom.domain.SpecificStorage.Value 1.0e-6

#-----------------------------------------------------------------------------
# Phases
#-----------------------------------------------------------------------------

pfset Phase.Names "water"

pfset Phase.water.Density.Type	Constant
pfset Phase.water.Density.Value	1.0

pfset Phase.water.Viscosity.Type	Constant
pfset Phase.water.Viscosity.Value	1.0

#-----------------------------------------------------------------------------
# Contaminants
#-----------------------------------------------------------------------------
pfset Contaminants.Names			""


#-----------------------------------------------------------------------------
# Gravity
#-----------------------------------------------------------------------------

pfset Gravity				1.0

#-----------------------------------------------------------------------------
# Setup timing info
#-----------------------------------------------------------------------------
 
pfset TimingInfo.BaseUnit        1.0
pfset TimingInfo.StartCount      0
pfset TimingInfo.StartTime       0.0
pfset TimingInfo.StopTime        5
pfset TimingInfo.DumpInterval    -1
pfset TimeStep.Type              Constant
pfset TimeStep.Value             1.0
 

#-----------------------------------------------------------------------------
# Porosity
#-----------------------------------------------------------------------------

pfset Geom.Porosity.GeomNames          domain

pfset Geom.domain.Porosity.Type    Constant
pfset Geom.domain.Porosity.Value   0.390

#-----------------------------------------------------------------------------
# Domain
#-----------------------------------------------------------------------------
pfset Domain.GeomName domain

#-----------------------------------------------------------------------------
# Mobility
#-----------------------------------------------------------------------------
pfset Phase.water.Mobility.Type        Constant
pfset Phase.water.Mobility.Value       1.0

#-----------------------------------------------------------------------------
# Relative Permeability
#-----------------------------------------------------------------------------
 
pfset Phase.RelPerm.Type               VanGenuchten
pfset Phase.RelPerm.GeomNames          "domain"
 
pfset Geom.domain.RelPerm.Alpha         3.5
pfset Geom.domain.RelPerm.N             2.

#---------------------------------------------------------
# Saturation
#---------------------------------------------------------

pfset Phase.Saturation.Type              VanGenuchten 
pfset Phase.Saturation.GeomNames         "domain"
 
pfset Geom.domain.Saturation.Alpha        3.5
pfset Geom.domain.Saturation.N            2.
pfset Geom.domain.Saturation.SRes         0.01
pfset Geom.domain.Saturation.SSat         1.0

#-----------------------------------------------------------------------------
# Wells
#-----------------------------------------------------------------------------
pfset Wells.Names ""


#-----------------------------------------------------------------------------
# Time Cycles
#-----------------------------------------------------------------------------
pfset Cycle.Names constant
pfset Cycle.constant.Names		"alltime"
pfset Cycle.constant.alltime.Length	 1
pfset Cycle.constant.Repeat		-1

#-----------------------------------------------------------------------------
# Boundary Conditions: Pressure
#-----------------------------------------------------------------------------
pfset BCPressure.PatchNames                   [pfget Geom.domain.Patches]
 
pfset Patch.x-lower.BCPressure.Type                   FluxConst
pfset Patch.x-lower.BCPressure.Cycle                  "constant"
pfset Patch.x-lower.BCPressure.alltime.Value          0.0
 
pfset Patch.y-lower.BCPressure.Type                   FluxConst
pfset Patch.y-lower.BCPressure.Cycle                  "constant"
pfset Patch.y-lower.BCPressure.alltime.Value          0.0
 
pfset Patch.z-lower.BCPressure.Type                   FluxConst
pfset Patch.z-lower.BCPressure.Cycle                  "constant"
pfset Patch.z-lower.BCPressure.alltime.Value          0.0
 
pfset Patch.x-upper.BCPressure.Type                   FluxConst
pfset Patch.x-upper.BCPressure.Cycle                  "constant"
pfset Patch.x-upper.BCPressure.alltime.Value          0.0
 
pfset Patch.y-upper.BCPressure.Type                   FluxConst
pfset Patch.y-upper.BCPressure.Cycle                  "constant"
pfset Patch.y-upper.BCPressure.alltime.Value          0.0
 
pfset Patch.z-upper.BCPressure.Type                   OverlandFlow
##pfset Patch.z-upper.BCPressure.Type                FluxConst 
pfset Patch.z-upper.BCPressure.Cycle                  "constant"
pfset Patch.z-upper.BCPressure.alltime.Value          0.0

#---------------------------------------------------------
# Topo slopes in x-direction
#---------------------------------------------------------
 
pfset TopoSlopesX.Type "Constant"
pfset TopoSlopesX.GeomNames "domain"
pfset TopoSlopesX.Geom.domain.Value -0.001
 
#---------------------------------------------------------
# Topo slopes in y-direction
#---------------------------------------------------------
 
pfset TopoSlopesY.Type "Constant"
pfset TopoSlopesY.GeomNames "domain"
pfset TopoSlopesY.Geom.domain.Value 0.001
 
#---------------------------------------------------------
# Mannings coefficient 
#---------------------------------------------------------
 
pfset Mannings.Type "Constant"
pfset Mannings.GeomNames "domain"
pfset Mannings.Geom.domain.Value 5.52e-6

#-----------------------------------------------------------------------------
# Phase sources:
#-----------------------------------------------------------------------------

pfset PhaseSources.water.Type                         Constant
pfset PhaseSources.water.GeomNames                    domain
pfset PhaseSources.water.Geom.domain.Value        0.0
 
#-----------------------------------------------------------------------------
# Exact solution specification for error calculations
#-----------------------------------------------------------------------------
 
pfset KnownSolution                                      NoKnownSolution

#-----------------------------------------------------------------------------
# Set solver parameters
#-----------------------------------------------------------------------------
 
pfset Solver                                             Richards
pfset Solver.MaxIter                                     500
 
pfset Solver.Nonlinear.MaxIter                           15
pfset Solver.Nonlinear.ResidualTol                       1e-9
pfset Solver.Nonlinear.EtaChoice                         EtaConstant
pfset Solver.Nonlinear.EtaValue                          0.01
pfset Solver.Nonlinear.UseJacobian                       False
pfset Solver.Nonlinear.DerivativeEpsilon                 1e-12
pfset Solver.Nonlinear.StepTol                           1e-20
pfset Solver.Nonlinear.Globalization                     LineSearch
pfset Solver.Linear.KrylovDimension                      15
pfset Solver.Linear.MaxRestart                           2
 
pfset Solver.Linear.Preconditioner                       PFMGOctree

pfset Solver.Linear.Preconditioner.PFMGOctree.BoxSizePowerOf2 2
 
pfset Solver.LSM                                         CLM
pfset Solver.WriteSiloCLM                                True
pfset Solver.CLM.MetForcing                              1D
pfset Solver.CLM.MetFileName                             narr_1hr.sc3.txt.0
pfset Solver.CLM.MetFilePath                             ./


pfset Solver.WriteSiloEvapTrans                          True
pfset Solver.WriteSiloOverlandBCFlux                     True


# Initial conditions: water pressure
#---------------------------------------------------------
 
pfset ICPressure.Type                                   HydroStaticPatch
pfset ICPressure.GeomNames                              domain
pfset Geom.domain.ICPressure.Value                      -2.0
 
pfset Geom.domain.ICPressure.RefGeom                    domain
pfset Geom.domain.ICPressure.RefPatch                   z-upper



set num_processors [expr [pfget Process.Topology.P] * [pfget Process.Topology.Q] * [pfget Process.Topology.R]]
for {set i 0} { $i <= $num_processors } {incr i} {
    file delete drv_vegm.dat.$i
    file copy  drv_vegm.dat drv_vegm.dat.$i
    file delete drv_clmin.dat.$i
    file copy drv_clmin.dat drv_clmin.dat.$i
}

#-----------------------------------------------------------------------------
# Run and Unload the ParFlow output files
#-----------------------------------------------------------------------------


pfrun clm 
pfundist clm 

#
# Tests 
#
source ../pftest.tcl
set passed 1

if ![pftestFile clm.out.perm_x.pfb "Max difference in perm_x" $sig_digits] {
    set passed 0
}
if ![pftestFile clm.out.perm_y.pfb "Max difference in perm_y" $sig_digits] {
    set passed 0
}
if ![pftestFile clm.out.perm_z.pfb "Max difference in perm_z" $sig_digits] {
    set passed 0
}

for {set i 0} { $i <= 5 } {incr i} {
    set i_string [format "%05d" $i]
    if ![pftestFile clm.out.press.$i_string.pfb "Max difference in Pressure for timestep $i_string" $sig_digits] {
    set passed 0
    }
    if ![pftestFile clm.out.satur.$i_string.pfb "Max difference in Saturation for timestep $i_string" $sig_digits] {
    set passed 0
    }
}

set mask [pfload clm.out.mask.pfb]
set top [Parflow::pfcomputetop $mask]

pfsave $top -pfb "clm.out.top_index.pfb"

set data [pfload clm.out.press.00000.pfb]
set top_data [Parflow::pfextracttop $top $data]

pfsave $data -pfb "clm.out.press.00000.pfb"
pfsave $top_data -pfb "clm.out.top.press.00000.pfb"

pfdelete $mask
pfdelete $top
pfdelete $data
pfdelete $top_data

if ![pftestFile clm.out.top_index.pfb "Max difference in top_index" $sig_digits] {
    set passed 0
}

if ![pftestFile clm.out.top.press.00000.pfb "Max difference in top_clm.out.press.00000.pfb" $sig_digits] {
    set passed 0
}



if $passed {
    puts "clm : PASSED"
} {
    puts "clm : FAILED"
}

