(C++ Simulation Code)
#include<iostream.h>
#include<iomanip.h>
#include<fstream.h>
#include<time.h>
#include<stdlib.h>
#include<math.h>
const int x_grid=225,y_grid=100;
double f_0 = (-1.34*(.0001)),
//coriolis parameter at reference latitude
beta =
(4.49*(.0000000001)), //beta
plane approximation
gravity = .0249,
//gravity on Jupiter
length = 90000.,
//length of GRS
height = 40000.,
//height of GRS
radius = 71400.,
//radius of jupiter
diffusion=.0001,
//diffusion term
h_init =30.0,
//initial thickness of GRS
ch_x = length/x_grid,
//physical distance of each x-grid
space
ch_y = height/y_grid,
//physical distance of each y-grid
space
strat =
(2200.*2200.*f_0*f_0)/(gravity*h_init), //stratification term
pi;
//arrays that hold the velocities
and thickness for three different time levels
double u_old
[x_grid+1][y_grid+2],u_curr[x_grid+1][y_grid+2];
double u_new [x_grid+1][y_grid+2],v_old [x_grid+1][y_grid+2];
double v_curr[x_grid+1][y_grid+2],v_new [x_grid+1][y_grid+2];
double h_old [x_grid+1][y_grid+2],h_curr[x_grid+1][y_grid+2];
double h_new [x_grid+1][y_grid+2];
double coriolis[y_grid+2],
//array holding coriolis force at different latitudes
u_2[y_grid+2];
//array holding background velocities at different latitudes
/*The initial velocities and
thicknesses are put in the arrays. The Gaussian
hill(valley) is used for the thicknesses which get the vortex to
form. */
void init_array()
{
double xk,
yk,
rpolar,
fr,
y;
//distance
from the equation
int i,j,ip1,im1; /*ip1 is i plus one; used to reference term to
right of the
current term. im1 is i minus one; used to reference term
to left of thecurrent term.*/
pi = 4.*atan(1.0);
for (j=1;j<=y_grid;j++)
{
y = -0.5*height + ((j-1.)*ch_y);
//calculates distance from equator
u_2[j] =
.07*sin((j-1.)/(y_grid-1.)*2.*pi); /*calculate u_2 and place's it
in the apporpriate place in
the
array*/
coriolis[j] = 2.0*(f_0+beta*y);
//the appropriate space in the array
}
for (i=1;i<=x_grid;i++)
{
xk=(i-1.)/(x_grid-1.)-.5;
//calculate xk
for (j=1;j<=y_grid;j++)
{
yk =
(j-1.)/(y_grid-1.)-.5; //calculate
yk
rpolar =
sqrt(xk*xk+yk*yk)/.5;
//calculate rpolar
fr =
.01*exp(-(rpolar/.1)*(rpolar/.1)); //calculate fr
h_curr[i][j] =
h_init*(1.+fr); /*our Gaussian hill
that gives
us the initial form of the vortex. It
give's the thicknesses current time step
an initial value that is slightly
different than h_init*/
}
h_curr[i][0] = h_curr[i][2]; //all thicknesses
at j=0 equal to j=2
h_curr[i][y_grid+1] =
h_curr[i][y_grid-1]; /*all
thicknesses at j=y_grid+1
equal to thicknesses at j=y_grid-1*/
}
for (i=1;i<=x_grid;i++)
{
for (j=1;j<=y_grid;j++)
{
/*if i is at an
endpoint ip1 and im1 is going to equal the other endpoint*/
if (i==1) im1=x_grid;
else im1=i-1;
if (i==x_grid) ip1=1;
else ip1=i+1;
/* Makes the current east/west velocities equal
to the background velocity
with a small pertubation dependent on latitude*/
u_curr[i][j] = u_2[j]
- (gravity*strat/coriolis[j])*
(h_curr[i][j+1]-h_curr[i][j-1])/ch_y;
/* Makes the current north/south velocities
equal to a small number
dependent on
latitude*/
v_curr[i][j] =
(gravity*strat/coriolis[j])*
(h_curr[ip1][j]-h_curr[im1][j])/ch_x;
}
}
// This loop handles the boundary
conditions
for (i=1;i<=x_grid;i++)
{
u_curr[i][0] = u_curr[i][2]; /*makes all north/south velocities
at j=0 equal j=2*/
u_curr[i][y_grid+1] = u_curr[i][y_grid-1]; /* makes all north/south
velocities at j=y_grid+1 equal to
the velocities at
j=y_grid-1 */
v_curr[i][1] = 0; //makes
all north/south velocities at j=1 equal to 0
v_curr[i][y_grid] = 0; //makes all v velocities at j=y_grid equal to 0
v_curr[i][0] = -v_curr[i][2]; /*makes all v velocities at j=0
equal to the negative velocities at
j=2*/
v_curr[i][y_grid+1] =
-v_curr[i][y_grid-1]; /* makes
all v velocities at
j=y_grid+1 equal to the negative
velocities at
j=y_grid-1*/
/*Makes the
current thicknesses at j=y_grid+1 equal thicknesses at
j=y_grid-1 minus
a small perturbation. The perturbation is used to let the
thickness change with the flow but still be
constrained.*/
h_curr[i][y_grid+1] = h_curr[i][y_grid-1] -
ch_y*coriolis[y_grid]/(gravity*strat)*(u_curr[i][y_grid]-u_2[y_grid]);
/* Makes the
current thicknesses at j=0 equal to the thicknesses at j=2
minus a small perturbation.*/
h_curr[i][0] = h_curr[i][2] +
ch_y*coriolis[1]/(gravity*strat)*
(u_curr[i][1]-u_2[1]);
}
//Places the current time level
arrays into the older time level array
for (i=1;i<=x_grid;i++)
{
for (j=0;j<=y_grid+1;j++)
{
u_old[i][j]=u_curr[i][j];
v_old[i][j]=v_curr[i][j];
h_old[i][j]=h_curr[i][j];
}
}
}
/* This nested loop calculates the
new u velocity and stores it to u_new*/
void u_change(double ch_time,int ntime)
{
int i,j,ip1,im1;
//the equation was split into
elements so as not to confuse the parenthesis
double element1,element2,element3,element4,element5,element6;
for (i=1;i<=x_grid;i++)
{
for (j=1;j<=y_grid;j++)
{
if (i==1) im1=x_grid;
else im1=i-1;
if (i==x_grid) ip1=1;
else ip1=i+1;
element1=-u_curr[i][j]*((u_curr[ip1][j]-u_curr[im1][j])/ch_x);
element2=-v_curr[i][j]*((u_curr[i][j+1]-u_curr[i][j-1])/ch_y);
element3=coriolis[j]*v_curr[i][j];
element4=-gravity*strat*((h_curr[ip1][j]-h_curr[im1][j])/ch_x);
element5=(u_curr[ip1][j]-2.*u_curr[i][j]+u_curr[im1][j])/(ch_x*ch_x);
element6=(u_curr[i][j+1]-2.*u_curr[i][j]+u_curr[i][j-1])/(ch_y*ch_y);
if
(ntime==1) /* the
first time through ch_time has to be divided by 2
because the first time though we aren't making a full
leap frog because we don't have the new time level yet.*/
{
u_new[i][j]=(element1+element2+element3+element4+
diffusion*(element5+element6))*(ch_time/2.) +
u_old[i][j];
}
else
{
u_new[i][j]=(element1+element2+element3+element4+
diffusion*(element5+element6))*(ch_time) +
u_old[i][j]; //see finite difference
equation sheet for u
}
}
}
//boundary
conditions for u
for (i=1;i<=x_grid;i++)
{
u_new[i][0]=u_new[i][2]; //makes the east-west velocity at the next time
//level at j=0 equal to the
east-west velocity
//at j=2
u_new[i][y_grid+1]=u_new[i][y_grid-1];
//makes the
east-west velocity at the next time level at the bottom equal
} //to the east-west velocity two
rows up from the bottom
}
/* This loop calculates the new v
velocity and stores it to v_new*/
void v_change(double ch_time,int ntime)
{
double element1,element2,element3,element4,element5,element6;
int jp1,jm1,ip1,im1,i,j;
for (i=1;i<=x_grid;i++)
{
for (j=1;j<=y_grid;j++)
{
if (i==1) im1=x_grid;
else im1=i-1;
if (i==x_grid) ip1=1;
else ip1=i+1;
jm1=j-1;
jp1=j+1;
element1=-u_curr[i][j]*((v_curr[ip1][j]-v_curr[im1][j])/ch_x);
element2=-v_curr[i][j]*((v_curr[i][jp1]-v_curr[i][jm1])/ch_y);
element3=-coriolis[j]*(u_curr[i][j]-u_2[j]);
element4=-gravity*strat*((h_curr[i][jp1]-h_curr[i][jm1])/ch_y);
element5=(v_curr[ip1][j]-2.*v_curr[i][j]+v_curr[im1][j])/(ch_x*ch_x);
element6=(v_curr[i][jp1]-2.*v_curr[i][j]+v_curr[i][jm1])/(ch_y*ch_y);
if (ntime==1)
{
v_new[i][j]=(element1+element2+element3+element4+
diffusion*(element5+element6))*(ch_time/2.)+
v_old[i][j]; //see corresponding
section in u_change
}
else
{
v_new[i][j]=(element1+element2+element3+element4+
diffusion*(element5+element6))*(ch_time)+
v_old[i][j];
}
}
}
for(i=1;i<=x_grid;i++)
{
v_new[i][1]=0;
v_new[i][0]=-v_new[i][2];
v_new[i][y_grid]=0;
v_new[i][y_grid+1]=-v_new[i][y_grid-1];
}
}
void h_change(double ch_time,int ntime)
{
int i,j,im1,ip1,jp1,jm1;
double element1,element2,element3,element4;
for (i=1;i<=x_grid;i++)
{
for (j=1;j<=y_grid;j++)
{
if (i==1) im1=x_grid;
else im1=i-1;
if (i==x_grid) ip1=1;
else ip1=i+1;
jm1=j-1;
jp1=j+1;
element1=-(h_curr[ip1][j]*u_curr[ip1][j]-
h_curr[im1][j]*u_curr[im1][j])/ch_x;
element2=-(h_curr[i][jp1]*v_curr[i][jp1]-
h_curr[i][jm1]*v_curr[i][jm1])/ch_y;
element3=(h_curr[ip1][j]-2*h_curr[i][j]+h_curr[im1][j])/(ch_x*ch_x);
element4=(h_curr[i][jp1]-2*h_curr[i][j]+h_curr[i][jm1])/(ch_y*ch_y);
if (ntime==1)
{
h_new[i][j]=(element1+element2+diffusion*(element3+element4))*
(ch_time/2.)+h_old[i][j];
//see corresponding section in u_change
}
else
{
h_new[i][j]=(element1+element2+diffusion*(element3+element4))*
ch_time + h_old[i][j];
}
}
}
for(i=1;i<=x_grid;i++)
{
h_new[i][y_grid+1] = h_new[i][y_grid-1] -
(ch_y*coriolis[y_grid]/(gravity*strat)*
(u_new[i][y_grid]-u_2[y_grid]));
/*Makes the new
thicknesses at j=y_grid+1 equal thicknesses at
j=y_grid-1 minus a small perturbation.*/
h_new[i][0] =
h_new[i][2]+(ch_y*coriolis[1]/(gravity*strat)*
(u_new[i][1]-u_2[1]));
/*Makes the new
thicknesses at j= equal thicknesses at
j=2 minus a small perturbation.*/
}
}
void updat_array() //advances the
program one time level by making
{
//new values current,
and current values old.
int i,j;
for (i=1;i<=x_grid;i++)
{
for (j=0;j<=y_grid+1;j++)
{
u_old[i][j]=u_curr[i][j]; //makes
the current values for u, h, and v
v_old[i][j]=v_curr[i][j]; //become
the old values, to progress the
h_old[i][j]=h_curr[i][j]; //
program in time.
}
}
for (i=1;i<=x_grid;i++)
{
for (j=0;j<=y_grid+1;j++)
{
u_curr[i][j]=u_new[i][j]; //makes
the new values for u, h, and v become
v_curr[i][j]=v_new[i][j]; //the
current values for the next loop through
h_curr[i][j]=h_new[i][j]; //the
program
}
}
}
void main()
{
ofstream fout;
//set the precision to six places
after the decimal
fout.setf(ios::fixed,ios::floatfield);
fout.setf(ios::showpoint);
fout<<setprecision(6);
double ch_time; //change in the time
for each time through the loop
double temp, //used as temporary variable to store velocities
vmax; /*When we take the absolute value of the
velocities we store the largest number to
vmax.*/
int ntime,i,j;
init_array();
vmax = -1.0; //initialize vmax
for (i=1;i<=x_grid;i++)
{
for (j=1;j<=y_grid;j++)
{
if (u_curr[i][j]<0) temp=-u_curr[i][j];
else temp= u_curr[i][j];
if (temp>vmax) vmax=temp;
}
}
//ch_time is
calculated two different ways and the smallest is chosen
if (.1*((length/x_grid)/vmax)<(.1*(height/y_grid)/vmax))
ch_time=0.1*((length/x_grid)/vmax);
else ch_time=0.1*((height/y_grid)/vmax);
for (ntime=1;ntime<=200000;ntime++)
{
u_change(ch_time,ntime); //function
to calculate new u velocity
v_change(ch_time,ntime); //function
to calculate new v velocity
h_change(ch_time,ntime); //function
to calculate new thicknesses
updat_array();
//function to move the arrays up a
time level
/*in the following code certain time
steps is saved to data files. In some runs of the code the data
files are saved to datakdl and datamlf files because we changed
the Gaussian hill(valley). This code saves it to datambm files */
if (ntime==1)
{ fout.open("datambm1.dat"); //open file for initial values storage
for (j=1;j<=y_grid;j++)
{
for (i=1;i<=x_grid;i++)
{
fout<<" "<<v_new[i][j]<<"
"<<u_new[i][j]-u_2[j] //output data
<<" "<<h_new[i][j]<<"
"<<u_2[j]<<endl;
}
}
fout.close(); //close file
}
if (ntime==5000)
{ fout.open("datambm5k.dat"); //open file for 5,000 time steps
for (j=1;j<=y_grid;j++)
{
for (i=1;i<=x_grid;i++)
{
fout<<" "<<v_new[i][j]<<"
"<<u_new[i][j]-u_2[j] //output data
<<" "<<h_new[i][j]<<endl;
}
}
fout.close(); //close file
}
if (ntime==20000)
{ fout.open("datambm20k.dat"); //open file for 20,000 time steps
for (j=1;j<=y_grid;j++)
{
for (i=1;i<=x_grid;i++)
{
fout<<" "<<v_new[i][j]<<"
"<<u_new[i][j]-u_2[j] //output data
<<" "<<h_new[i][j]<<endl;
}
}
fout.close(); //close file
}
if (ntime==40000)
{ fout.open("datambm40k.dat"); //open file for 40,000 time steps
for (j=1;j<=y_grid;j++)
{
for (i=1;i<=x_grid;i++)
{
fout<<" "<<v_new[i][j]<<"
"<<u_new[i][j]-u_2[j] //output
data
<<" "<<h_new[i][j]<<endl;
}
}
fout.close(); //close file
}
if (ntime==60000)
{ fout.open("datambm60k.dat"); //open file for 60,000 time steps
for (j=1;j<=y_grid;j++)
{
for (i=1;i<=x_grid;i++)
{
fout<<" "<<v_new[i][j]<<"
"<<u_new[i][j]-u_2[j] //output
data
<<" "<<h_new[i][j]<<endl;
}
}
fout.close(); //close file
}
if (ntime==80000)
{ fout.open("datambm80k.dat"); //open file for 80,000 time steps
for (j=1;j<=y_grid;j++)
{
for (i=1;i<=x_grid;i++)
{
fout<<" "<<v_new[i][j]<<"
"<<u_new[i][j]-u_2[j] //output
data
<<" "<<h_new[i][j]<<endl;
}
}
fout.close(); //close file
}
if (ntime==100000)
{ fout.open("datambm100k.dat"); //open file for 100,000 time steps
for (j=1;j<=y_grid;j++)
{
for (i=1;i<=x_grid;i++)
{
fout<<" "<<v_new[i][j]<<"
"<<u_new[i][j]-u_2[j] //output data
<<" "<<h_new[i][j]<<endl;
}
}
fout.close(); //close file
}
if (ntime==120000)
{ fout.open("datambm120k.dat"); //open file for 120,000 time steps
for (j=1;j<=y_grid;j++)
{
for (i=1;i<=x_grid;i++)
{
fout<<" "<<v_new[i][j]<<"
"<<u_new[i][j]-u_2[j] //output
data
<<" "<<h_new[i][j]<<endl;
}
}
fout.close(); //close file
}
if (ntime==140000)
{ fout.open("datambm140k.dat"); //open file for 140,000 time steps
for (j=1;j<=y_grid;j++)
{
for (i=1;i<=x_grid;i++)
{
fout<<" "<<v_new[i][j]<<"
"<<u_new[i][j]-u_2[j] //output
data
<<" "<<h_new[i][j]<<endl;
}
}
fout.close(); //close file
}
if (ntime==160000)
{ fout.open("datambm160k.dat"); //open file for 160,000 time steps
for (j=1;j<=y_grid;j++)
{
for (i=1;i<=x_grid;i++)
{
fout<<" "<<v_new[i][j]<<"
"<<u_new[i][j]-u_2[j] //output
data
<<" "<<h_new[i][j]<<endl;
}
}
fout.close(); //close file
}
if (ntime==180000)
{ fout.open("datambm180k.dat"); //open file for 180,000 time steps
for (j=1;j<=y_grid;j++)
{
for (i=1;i<=x_grid;i++)
{
fout<<" "<<v_new[i][j]<<"
"<<u_new[i][j]-u_2[j] //output
data
<<" "<<h_new[i][j]<<endl;
}
}
fout.close(); //close file
}
if (ntime==200000)
{ fout.open("datambm200k.dat"); //open file for 200,000 time steps
for (j=1;j<=y_grid;j++)
{
for (i=1;i<=x_grid;i++)
{
fout<<" "<<v_new[i][j]<<"
"<<u_new[i][j]-u_2[j] //output data
<<" "<<h_new[i][j]<<endl;
}
}
fout.close(); //close file
}
}
}