Define LENGTH (150./180.*M_PI)

1 0 GfsRiver GfsBox GfsGEdge {} {
    PhysicalParams { L = LENGTH }
    MetricLonLat M 1.
    Refine 8
    InitFraction P (0.2 - acos(cos(x*M_PI/180.)*cos(y*M_PI/180.)))
    Init {} { P = 0.2 + 1.8*P/LENGTH }
    Time { end = 0.9 }
    OutputTime { istep = 10 } stderr
    GModule gfsview
    OutputView { start = 0.3 step = 0.3 } isolines-%g.eps { 
	format = EPS line_width = 0.5 
    } isolines.gfv
    OutputSimulation { step = 0.3 } sim-%g.txt { variables = U,V,P format = text }
    OutputScalarSum { istep = 1 } sp { v = P }
#    OutputSimulation { istep = 10 } stdout
    EventScript { start = end } {
	status=0
	for i in 0.3 0.6 0.9; do
	    if awk 'BEGIN { n1 = 0; pi = 3.14159265359 }{
              if ($1 != "#") {
                c = cos($1*pi/180.)*cos($2*pi/180.);
                d = atan2(sqrt(1. - c*c),c)*180./pi
                i = int(d*2.)
                x[i] += d
                y[i] += $6
                n[i]++
                x1[n1] = d;
                y1[n1++] = $6;
              }
            }END {
              for (i = 0; i <= 180; i++)
                if (n[i] > 0)
                  print x[i]/n[i], y[i]/n[i], n[i];
              sum = 0.;
              for (i = 0; i < n1; i++) {
                j = int(x1[i]*2.)
                d = y1[i] - y[j]/n[j];
                sum += d*d;
              }
              scatter = sqrt(sum/n1);
              if (scatter > 0.013) {
                print scatter > "/dev/stderr";
                exit (1);
              }
            }' < sim-$i.txt > prof-$i.txt ; then :
	    else
		status=$GFS_STOP;
	    fi

	    cat <<EOF | gnuplot
rdist(x,y)=acos(cos(x*pi/180.)*cos(y*pi/180.))*180./pi
cdist(x,y)=sqrt(x*x+y*y)
set xlabel 'Angular distance (degree)'
set ylabel 'Surface height'
set xtics 0,22.5,90
set ytics 0,0.25,0.75
set term postscript eps color lw 2 solid 20
set output 'p-$i.eps'
plot [0:90][0:0.75]'sim-$i.txt' u (rdist(\$1,\$2)):6 ps 0.25 pt 6 t '', 'prof-$i.txt' w l lw 2 t ''
EOF

	done
	exit $status
    }
}
GfsBox {
    right = Boundary { BcNeumann U 0 }
    left = Boundary { BcNeumann U 0 }
    top = Boundary { BcNeumann V 0 }
    bottom = Boundary { BcNeumann V 0 }
}