TALK=T;RUN(1,1)

  DISPLAY
  
  In-Form is used to supply, at the boundaries of a 2D box with x 
  horizontal and y vertical, boundary conditions of pressure and 
  velocity which correspond to an ideal-fluid travelling gravity 
  wave in moderately deep water.
  
  Within the box, the pressure and velocity are computed from the
  momentum and energy equations in the usual way.
  
  It is interesting to compare their values with those 
  corresponding to the ideal-fluid wave, namely: pdif, udif and vdif
  
  If OBSTCL=T, an obstacle appears on the the sea-bed.
  If STRA  =T, the stresses and strains in it are computed
   
  ENDDIS
   
  Constants used for connecting p1, u1, v1 with ppot, upot and vpot
  at the boundaries.
  Declarations:
real(pcow,ucow,vcow,pcoe,ucoe,vcoe)
integer(lstepw)  ! duration of whole-field patches
boolean(obstcl)
  Settings
STRA   = F
obstcl = f  
obstcl = t  
ucow=1.e10; vcow=1.e10; pcow=1.e10
ucoe=1.e10; vcoe=1.e10; pcoe=1.e5
 
lstepw=2 ! flow field is fixed to potential for first 2 time steps
  
nx=50      ! discretization in space
ny=50
  nx=5
  ny=5
lstep=50   ! and time
  lstep=5   ! and time

  lstep=100
             ! surface-wave-defining constants
real(lamda,aa,m,h,sig,g,arg,speed,const,twopi,tim1)
g=9.81       ! gravitational acceleration
h=40.0       ! total depth of water 
lamda=120    ! wave-length
  lamda=1200    ! wave-length
twopi=2.0*3.1416
m=twopi/lamda
arg=m*h
const=tanh(arg)
sig=(g*m*const)**0.5
speed=sig/m
tlast=lamda/speed
tlast=tlast*lstep/100

tim1=tlast/lstep ! time after first time interval
sig
speed
tlast
tim1
aa=3.0           ! amplitude of wave
  aa=0.3         ! amplitude of wave


xulast=lamda ! extent of domain in wave-travel direction
yvlast=0.75*h ! domain extent above ocean floor
steady=f   
resfac=0.
resref(1)=0.0
endit(1) =0.0
resref(3)=0.0
endit(3) =0.0
resref(5)=0.0
endit(5) =0.0
isg21=lsweep
#unigrid   ! use uniform grid
solve(p1,u1,v1)
rho1=1.000e3
  rho1=1.
  liter(p1)=1000
  store(imb1,pcor,min1)
  store(u1sl,v1sl,u1ad,v1ad)
(stored var udif is u1-upot )  
(stored var vdif is v1-vpot )  
(stored var pdif is p1-ppot ) 
store(pdcx) 
relax(u1,falsdt,1.e6)
relax(v1,falsdt,1.e6)

  !!!! patches !!!!!
  !!!! for inform-created variables !!!!!
patch(whole,phasem,1,nx,1,ny,1,nz,1,lstep)   ! for pot, ppot, upot, vpot
                                             ! and initial values
  !!!! zero-pressure outer patches
patch(p1x,phasem,1,1,1,ny,1,nz,1,lstep)      ! west edge
coval(p1x,p1,fixval,0.0)

patch(pny,phasem,1,nx,ny,ny,1,nz,1,lstep)    ! north edge
coval(pny,p1,fixval,0.0)

patch(pnx,phasem,nx,nx,1,ny,1,nz,1,lstep)    ! east edge
coval(pnx,p1,fixval,0.0)

  !!!! next-inner pressure patches
patch(>ppot2x,phasem,2,2,1,ny-1,1,nz,1,lstep)        ! near west edge
coval(>ppot2x,p1,pcoe,1)  

patch(>ppotny1,phasem,2,nx-1,ny-1,ny-1,1,nz,1,lstep) ! near north edge
coval(>ppotny1,p1,pcoe,1)  

patch(>ppotnx1,phasem,nx-1,nx-1,1,ny-1,1,nz,1,lstep)   ! near east edge
coval(>ppotnx1,p1,pcoe,1)  
  

  !!!! velocity patches near boundaries
patch(>upot1x,phasem,1,1,1,ny,1,nz,1,lstep)          ! near west edge
coval(>upot1x,u1,ucoe,1)  

patch(>vpot1x,phasem,1,1,1,ny-1,1,nz,1,lstep)  
coval(>vpot1x,v1,vcoe,1)  

patch(>upotny,phasem,1,nx,ny,ny,1,nz,1,lstep)        ! near north edge
coval(>upotny,u1,ucoe,1)      

patch(>vpotny1,phasem,1,nx,ny-1,ny-1,1,nz,1,lstep)
coval(>vpotny1,v1,vcoe,1)       ! edges       

patch(>upotnx1,phasem,nx-1,nx-1,1,ny-1,1,nz,1,lstep) ! near east edge
coval(>upotnx1,u1,ucoe,1)  

patch(>vpotnx,phasem,nx,nx,1,ny-1,1,nz,1,lstep)
coval(>vpotnx,v1,vcoe,1)  
  
  !!!! whole-field patches for lstepw time steps
patch(>ppotw,phasem,1,nx,1,ny,1,nz,1,lstepw)   
coval(>ppotw,p1,pcow ,1)  

patch(>upotw,phasem,1,nx,1,ny,1,nz,1,lstepw)   
coval(>upotw,u1,ucow,1)        ! whole domain

patch(>vpotw,phasem,1,nx,1,ny-1,1,nz,1,lstepw)  
coval(>vpotw,v1,vcow,1)        ! whole domain

               ! defining the velocity potential and derived quantities
char(form)     ! declare character variable
               
  
  Formula for the potential as function of space and time
form=aa*(cosh(m*yg))*cos(m*xg-sig*tim)
(stored var pot is :form:) 
form=aa*(cosh(m*yg))*cos(m*xg-sig*tim1)  ! note use of tim1
(initial of pot at whole is :form:)
  
  Formula for the potential-derived u velocity =
  - d pot/dx     Note: d/dx of cos mx = - m*sin ; so plus is OK
form = aa*m*(cosh(m*yg))*sin(m*xu-sig*tim) ! for all times
(stored var upot at whole is :form:)
form = aa*m*(cosh(m*yg))*sin(m*xu-sig*tim1) ! note use of tim1 here
(initial of upot at whole is :form:)
(initial of u1 at whole is :form:)

  Formula for the potential-derived v velocity, =
   - d pot / dy  Note: d/dy of cosh my is sinh my ; so minus is appropriate
form = - aa*m*(sinh(m*yv))*cos(m*xg-sig*tim)
(stored var vpot at whole is :form:)
form = - aa*m*(sinh(m*yv))*cos(m*xg-sig*tim1)  ! note use of tim1
(initial of vpot at whole is :form:)
(initial of v1 at whole is :form:)

  Formula for potential-derived  pressure =
  d pot / dt     Note: d/dt of cos mx - sig t = -(-sig)*sin ; so plus is OK 
const=sig*aa*rho1   
form = const*(cosh(m*yg))*sin(m*xg-sig*tim)
(stored var ppot at whole is :form:)
form = const*(cosh(m*yg))*sin(m*xg-sig*tim1) ! note use of tim1
(initial of ppot at whole is :form:)
(initial of p1 at whole is :form:)
if(nz.eq.1) then
idispa=lstep/50       ! for dumping to parphi
endif
text(p :pcow: :pcoe: u :ucow: :ucoe:

    obstacle to u1 on floor
  patch(obstcl1,phasem,nx/2,nx/2,1,ny/2,1,1,1,lstep)
  patch(obstcl1,phasem,nx/2,nx/2,1,ny/2,1,1,1,lstep)
  coval(obstcl1,u1,1.e10,0.0)
                          
                          
                          
                          
                  structure on ocean floor        
                         __________
                C        !  ____  !
                B        !  !//!  ! 
                         !  !//!  ! 
                         !  !//!  ! 
                         !  !//!  ! 
                A        !  !//!  ! 
                            !//!  ! 
  --------------------------!//!--!-----------------------------------
                         1  2  3  4
                         
integer(ix1,ix2,ix3,ix4,iyA,iyB,iyC)
ix1=nx/2
ix2=ix1+2
ix3=ix2
ix4=ix3+2
iyA=5
iyB=iyA+10
iyC=iyB+5
                         
store(prps)
fiinit(prps)= -1
if(obstcl) then
iya=1
patch(wall11,inival,ix1,ix1,iyA,iyC,1,1,1,lstep)
init(wall11,prps,0,100)

patch(wall23,inival,ix2,ix3,1,iyB,1,1,1,lstep)
init(wall23,prps,0,100)

patch(wall44,inival,ix4,ix4,1,iyC,1,1,1,lstep)
init(wall44,prps,0,100)

patch(topC,inival,ix1,ix4,iyC+1,iyC+1,1,1,1,lstep)
init(topC,prps,0,100)

patch(bottom,phasem,ix3,ix4,1,1,1,1,1,lstep)
coval(bottom,p1,fixval,0.0)
endif
#maxabs                   
                          
  # the number of spedat special-print commands
SPEDAT(PRINT,NUMBER,I,3)  
                          
  # the command to print minimum and maximum values 
SPEDAT(PRINT,COMMAND1,C,MINMAX_p1)
SPEDAT(PRINT,COMMAND2,C,MINMAX_u1)
SPEDAT(PRINT,COMMAND3,C,MINMAX_v1)
  SPEDAT(PRINT,COMMAND4,C,MINMAX_v1sl)
  SPEDAT(PRINT,COMMAND5,C,MINMAX_v1ad)
  SPEDAT(PRINT,COMMAND6,C,MINMAX_imb1)
  SPEDAT(PRINT,COMMAND7,C,MINMAX_pcor)
   -----------------------------------------
 if(stra) then
  Additions for stress
 STORE(EPSZ,EPSY,EPSX,STRZ,STRY,STRX,DRH1,DVO1,ENUL,DEN1,PRPS)
 PATCH(BOTTOM11,SOUTH ,   ix1,   ix1,   1,   1,   1,   1,   1,  50)
 COVAL(BOTTOM11,U1  , 1.000000E+10, 0.000000E+00)
 COVAL(BOTTOM11,V1  , 1.000000E+10, 0.000000E+00)
 
 PATCH(BOTTOM23,SOUTH ,   ix2,   ix3,   1,   1,   1,   1,   1,  50)
 COVAL(BOTTOM23,U1  , 1.000000E+10, 0.000000E+00)
 COVAL(BOTTOM23,V1  , 1.000000E+10, 0.000000E+00)
 
 PATCH(BOTTOM44,SOUTH ,   ix4,   ix4,   1,   1,   1,   1,   1,  50)
 COVAL(BOTTOM44,U1  , 1.000000E+10, 0.000000E+00)
 COVAL(BOTTOM44,V1  , 1.000000E+10, 0.000000E+00)
 SPEDAT(SET,RLXFAC,RLXU1D,R,7.50000E-01)
 SPEDAT(SET,LONGNAME,P1,C,PRESSURE_OR_DILATATION)
 SPEDAT(SET,LONGNAME,ENUL,C,VISCOSITY_OR_LPLUS2G*1.E11)
 SPEDAT(SET,LONGNAME,DRH1,C,COMPRESSIBILITY_OR_LPLUSG*1.E11)
 SPEDAT(SET,LONGNAME,DVO1,C,TH.EXP.CO_R_OR_H*1.E11)
 SPEDAT(SET,LONGNAME,U1,C,X-DIRECTION_VELOCITY_OR_DISPLACEM$)
 SPEDAT(SET,LONGNAME,U1,C,ENT)
 SPEDAT(SET,LONGNAME,EPSX,C,X-DIRECTION_STRAIN)
 SPEDAT(SET,LONGNAME,STRX,C,X-DIRECTION_STRESS)
 SPEDAT(SET,LONGNAME,V1,C,Y-DIRECTION_VELOCITY_OR_DISPLACEM$)
 SPEDAT(SET,LONGNAME,V1,C,ENT)
 SPEDAT(SET,LONGNAME,EPSY,C,Y-DIRECTION_STRAIN)
 SPEDAT(SET,LONGNAME,STRY,C,Y-DIRECTION_STRESS)
 SPEDAT(SET,LONGNAME,U1,C,X-DIRECTION_VELOCITY_OR_DISPLACEM$)
 SPEDAT(SET,LONGNAME,U1,C,ENT)
 SPEDAT(SET,LONGNAME,V1,C,Y-DIRECTION_VELOCITY_OR_DISPLACEM$)
 SPEDAT(SET,LONGNAME,V1,C,ENT)
 SPEDAT(SET,LONGNAME,EPSZ,C,Z-DIRECTION_STRAIN)
 SPEDAT(SET,LONGNAME,STRZ,C,Z-DIRECTION_STRESS)
 SPEDAT(SET,LONGNAME,EPST,C,LINEARTHERMALEXPANSION)
 SPEDAT(SET,PRINT,NUMBER,I,7)
 SPEDAT(SET,PRINT,COMMAND1,C,MINMAX_P1)
 SPEDAT(SET,PRINT,COMMAND2,C,MINMAX_U1)
 SPEDAT(SET,PRINT,COMMAND3,C,MINMAX_V1)
 SPEDAT(SET,PRINT,COMMAND3,C,MINMAX_strx)
 SPEDAT(SET,PRINT,COMMAND3,C,MINMAX_stry)
 SPEDAT(SET,PRINT,COMMAND3,C,MINMAX_epsx)
 SPEDAT(SET,PRINT,COMMAND3,C,MINMAX_epsy)
 endif                    ! end of stress additions
difcut=0
ixmon=nx/2;iymon=ny/2;izmon=nz/2
tstswp=-1
nxprin=1
nyprin=1
iyprf=1;iyprl=iyc
ixprf=ix1-2;ixprl=ix4-2
ntprin=lstep/5
  nprint=1
  
  informbegin
  debug t
  initial f
  store t
  formula t
  informend
  

  debug=t
lsweep=25
lsweep=10
  lsweep=1
  lstep=3
  ntprin=1
istprf=3
istdb1=3
istdb2=3
  iyprf=1
iyprl=2
  ixprf=1;ixprl=4
dbindx=t
  search=t;idbf0=108591
dbcomp=t
dbsoda=t
dbcmpn=t
dbcmpe=t
dbsol2=t
  dbgphi(1)=t
dbgphi(3)=t
  dbgphi(5)=t
  dbgphi(9)=t
  dbflux=t
  dbconv=t
stop