UPPER LEVEL FRONTOGENESIS

SetDirectory["D:/cygwin/home/legras/h2/cours2/M2-approf/fronto"]

<<ColorbarPlot`

Basic flow

References: Keyser & Pecnick, 1985, J. Atmos. Sci., 42, pp.1259-1282; Keyser & Pecnick, 1985, J. Atmos. Sci., 42, pp. 1283-1305; Keyser & Shapiro, 1986, Mon. Wea. Rev., 114, pp. 452-499; Hoskins & Bretherton, 1972, J. Atmos. Scio., 29, 11-37; Hoskins, 1982, Ann. Rev. Fluid Mech., 14, pp. 131-151.

Note: Solution exacte de l' équation du mouvement avec confluence et advection de température

u= - α x + γ(z-(Subscript[z, r])) (e^(α t))
v= α y
Φ = g z + α f x y - ((α^2)/2)((x^2)+(y^2)) -γ f y (z-(Subscript[z, r])) (e^(α t))
b = g - γ y f (e^(α t))

Basic state at initial time in the y=0 plane

h[x] : hauteur de la tropopause
b[x_, z_] : champ de température
Φ[x_, z_] : géopotentiel
Ug, Vg : vitesse géostrophique

h[x_]:=(Subscript[h, 0])+ K(((((Subscript[N, 2])^2)-((Subscript[N, 1])^2)))^(-1))ArcTan[(x/L)];
b[x_,ζ_]:= g + (((Subscript[N, 1])^2)ζ + K ArcTan[(x/L)])UnitStep[h[x]-ζ] + (((Subscript[N, 2])^2)ζ-(((Subscript[N, 2])^2)-((Subscript[N, 1])^2))(Subscript[h, 0]))UnitStep[ζ-h[x]];
Φ[x_,ζ_]:= g ζ  + K Arc.Tan[(x/L)] (ζ-(Subscript[z, c])) UnitStep[h[x]-ζ]  + ((((Subscript[N, 2])^2)-((Subscript[N, 1])^2)) ((1/2)(h^2)- h (Subscript[z, c]) - (Subscript[h, 0]) (ζ-(Subscript[z, c])))+ (1/2)((Subscript[N, 2])^2) (z^2))UnitStep[ζ-h[x]];
Ug[x_,ζ_]:=-α x + (γ/f)(ζ-(Subscript[z, r])) ;
Vg[x_,ζ_] := (K/(f L))(((1+((((x/L)))^2)))^(-1)) (ζ-(Subscript[z, c]))UnitStep[h[x]-ζ] + (K/(f L))(((1+((((x/L)))^2)))^(-1))(h[x]-(Subscript[z, c]))UnitStep[ζ-h[x]];
K= ((Δθ g)/(π (Subscript[θ, 0])));

Q1 component of the Q-vector and elements of the Sawyer-Eliassen equation

Q1[x_,ζ_]=Simplify[f ((∂_x)Vg[x,ζ] (∂_ζ)Ug[x,ζ] -(∂_x)Ug[x,ζ] (∂_ζ)Vg[x,ζ] )/.  DiracDelta[t_]->0];
Ns2[x_,ζ_]=Simplify[(∂_ζ)b[x,ζ] /.  DiracDelta[t_]->0];
S2[x_,ζ_]=Simplify[(∂_x)b[x,ζ]/.  DiracDelta[t_]->0];
F2[x_,ζ_] = Simplify[f(f+(∂_x)Vg[x,ζ]) /.  DiracDelta[t_]->0];
Ng2[x_,ζ_]=Simplify[((Subscript[N, 1])^2)UnitStep[h[x]-ζ] + ((Subscript[N, 2])^2) UnitStep[ζ-h[x]]];PV[x_,ζ_]=Simplify[ (1/f)(F2[x,ζ]Ns2[x,ζ] - ((S2[x,ζ])^2))];

Values of the parameters

Parameter list

f:	Coriolis parameter
g:	gravity
Lx, H: width and height of the computing domain
...

Pure confluence

distances are in km, times in second

param1={f->(10^(-4)),g->9.81 (10^(-3)),
(Subscript[h, 0])->8.15,H->13.5,
Lx->2000.,L->500.,
((Subscript[N, 1])^2) -> (10^(-4)),((Subscript[N, 2])^2) ->9. (10^(-4)),
(Subscript[z, r])->5.,(Subscript[z, c])->0.,(*after Keyser and Pecnick,1985*)
α->(10^(-5)),Δθ -> 40., (Subscript[θ, 0])->273,
γ->0 };
param=param1;

Confluence + cold advection

Keyser & Pecnick use L=796, try also L=525

param2={f->1. (10^(-4)),g->9.81 (10^(-3)),
(Subscript[h, 0])->8.15,H->13.5,
Lx->2000.,L->500.,
((Subscript[N, 1])^2) -> 1. (10^(-4)),((Subscript[N, 2])^2) ->9. (10^(-4)),
(Subscript[z, r])->5.,(Subscript[z, c])->0.,
α->1. (10^(-5)),Δθ -> 40., (Subscript[θ, 0])->273,
γ->-5.74  (10^(-7)) };
param=param2;

Confluence + warm advection

param3={f->1. (10^(-4)),g->9.81 (10^(-3)),
(Subscript[h, 0])->8.15,H->13.5,
Lx->2000.,L->500.,
((Subscript[N, 1])^2) -> 1. (10^(-4)),((Subscript[N, 2])^2) ->9. (10^(-4)),
(Subscript[z, r])->5.,(Subscript[z, c])->0.,
α->1. (10^(-5)),Δθ -> 40., (Subscript[θ, 0])->273,
γ->5.74  (10^(-7)) };
param=param3;

Plots

Horizontal temperature ang geostrophic wind [in the (x,y)-plane] and cross front temperature distribution [in the (x,z)-plane]

g1:=VectorPlot[Evaluate[{Ug[x,(Subscript[z, r])],Vg[x,(Subscript[z, r])]+α y}/.param],{x,-2000,2000},{y,0,4000},DisplayFunction->Identity];
g2:=ContourPlot[Evaluate[(((b[x,(Subscript[z, r])]-γ y) (Subscript[θ, 0]))/g)/.param],{x,-2000,2000},{y,0,4000.},Contours->50,ContourStyle->None,ColorFunction->(RGBColor[Min[2 ((#1)^2),1],4 #1 (1-#1),Min[2 (((1-#1))^2),1]]&),DisplayFunction->Identity]/.param;
g3:=ContourPlot[Evaluate[((b[x,z] (Subscript[θ, 0]))/g)/.param],{x,-1000,1000},{z,0,13.5},Contours->Table[-50+5 i,{i,200}],ColorFunction->Hue,DisplayFunction->Identity,FrameLabel->{"x(km)","z(km)"},Axes->True,PlotLabel->"Temperature",Exclusions->None];
g4:=Show[GraphicsRow[{Show[g2,g1,FrameLabel->{"x(km)","y(km)"},Axes->True,PlotLabel->"Wind + temperature"],g3}]]

g4

Along front geostrophic wind (x,z)-plane

ContourPlot[Evaluate[Vg[x,z]/.param],{x,-2000,2000},{z,0,13.5},Contours->10,ColorFunction->Hue,FrameLabel->{"x(km)","z(km)"},Axes->True,Exclusions->None]/.param

g1=ColorbarPlot[Evaluate[1000 Vg[#1,#2] /.param]&,{-2000,2000},{0,13.5},PlotType->"Contour",Axes->True,XLabel->"x (km)",YLabel->"z (km)",CLabel->"(m/s)",Contours->10,Title->"Vg along front",ColorFunction->"Rainbow" ,Exclusions->None,Height->400]

Cross front  geostorphic wind (x,z)-plane

g2=ContourPlot[Evaluate[Ug[x,z] /.param],{x,-2000,2000},{z,0,13.5},Contours->10,ColorFunction->Hue] /. param

g1=ColorbarPlot[Evaluate[1000 Ug[#1,#2] /.param]&,{-2000,2000},{0,13.5},PlotType->"Contour",Axes->True,XLabel->"x (km)",YLabel->"z (km)",CLabel->"(m/s)",Contours->10,Title->"Ug across front",ColorFunction->"Rainbow" ,Exclusions->None,Height->400]/. param

Cross front and along front geostorphic wind (x,z)-plane

g5:=Show[GraphicsRow[{ColorbarPlot[Evaluate[1000 Vg[#1,#2] /.param]&,{-1000,1000},{0,13.5},PlotType->"Contour",Axes->True,XLabel->"x (km)",YLabel->"z (km)",CLabel->"(m/s)",Contours->10,Title->"Vg along front",ColorFunction->"Rainbow" ,Exclusions->None,Height->400],
ColorbarPlot[Evaluate[1000 Ug[#1,#2] /.param]&,{-1000,1000},{0,13.5},PlotType->"Contour",Axes->True,XLabel->"x (km)",YLabel->"z (km)",CLabel->"(m/s)",Contours->10,Title->"Ug across front",ColorFunction->"Rainbow" ,Exclusions->None,Height->400]}]]

g5

Potential vorticity (x,z)-plane

g6:=Show[ContourPlot[Evaluate[PV[x,z]/.param],{x,-1000,1000},{z,0,13.5},ContourStyle->None,ColorFunction->Hue,DisplayFunction->Identity,Exclusions->None],ContourPlot[Evaluate[PV[x,z]/.param],{x,-1000,1000},{z,0,13.5},Contours->{0.},ContourShading->False,DisplayFunction->Identity,Exclusions->None],FrameLabel->{"x(km)","z(km)"},PlotLabel->"Potential vorticity",Axes->True]

g6

x-component of Q-vector (x,z)-plane

ContourPlot[Evaluate[Q1[x,z] /.param],{x,-2000,2000},{z,0,13.5},Contours->10,ColorFunction->"BrightBands",FrameLabel->{"x(km)","z(km)"},Axes->True,Exclusions->None] /. param

g7:=ColorbarPlot[Evaluate[Q1[#1,#2] /.param]&,{-1000,1000},{0,13.5},PlotType->"Contour",Axes->True,XLabel->"x (km)",YLabel->"z (km)",CLabel->"(1/s^2)",Contours->10,Title->"Q1",ColorFunction->"BrightBands" ,Exclusions->None,Height->400]/. param

g7

Solution of the QG Sawyer-Eliassen equation

SolveQGSE :=Module[{},num=100;hstep=2. Lx/(num+1)/. param; vstep= H/(num+1) /. param;
vars=Table[a[i,j],{j,num},{i,num}];
hkern2={1,-2,1}/hstep^2; 
vkern2={1,-2,1}/vstep^2; 
pQ1[i_,j_]=Simplify[N[Q1[-Lx + i hstep,j vstep] /. param]];
pNg2[i_,j_]=Simplify[N[Ng2[-Lx + i hstep, j vstep] /. param]];
mNg2 = Table[pNg2[i,j],{j,num},{i,num}];
mF2 = Table[N[(f^2)]/. param,{j,num},{i,num}];
lap=mNg2 Map[ListCorrelate[hkern2,#,{2,2},0]&,vars]+ mF2  Transpose[Map[ListCorrelate[vkern2,#,{2,2},0]&,Transpose[vars]]];
eqns=Thread[Map[Flatten,lap==Table[-2 pQ1[i,j],{j,num},{i,num}]]];
{vec,mat}=CoefficientArrays[eqns,Flatten[vars]];
sol2=Partition[LinearSolve[mat,-vec],num]; //Timing
]

MakePsi:= Module[{},ZeroFrame=Prepend[Append[#1,0],0]&;
sol=Transpose[ZeroFrame/@Transpose[ZeroFrame/@sol2]];
Ψ=ListInterpolation[Transpose[sol],{{-2000.,2000.},{0,1350.}}]]

g8:=ContourPlot[Ψ[x,z],{x,-1000.,1000.},{z,0,1350},ColorFunction->Hue,DisplayFunction->Identity];
g9:=VectorPlot[Evaluate[{1000 (∂_z)Ψ[x,z],- 1000 (∂_x)Ψ[x,z]}],{x,-1000.,1000.},{z,0.,1350},DisplayFunction->Identity];
g10:=Show[GraphicsRow[{Show[{g8,g9},FrameLabel->{"x(km)","z(dm)"}],g3}]]

SolveQGSE;
MakePsi; g10

Solution of the full Sawyer-Eliassen equation

SolveSE := Module[{},num=100;hstep=2. Lx/(num+1)/. param; vstep= H/(num+1) /. param;
vars=Table[a[i,j],{j,num},{i,num}];
hkern2={1,-2,1}/hstep^2; hkern1={1,0,-1}/(2 hstep);
vkern2={1,-2,1}/vstep^2; vkern1 ={1,0,-1}/(2 vstep);
pQ1[i_,j_]=Simplify[N[Q1[-Lx + i hstep,j vstep] /. param]];
pNg2[i_,j_]=Simplify[N[Ng2[-Lx + i hstep,j vstep] + Ns2[-Lx + i hstep,j vstep]/. param]];
pF2[i_,j_]=Simplify[N[F2[-Lx + i hstep,j vstep] /. param]];
pS2[i_,j_]=Simplify[N[S2[-Lx + i hstep,j vstep] /. param]];
mNg2 = Table[pNg2[i,j],{j,num},{i,num}];
mF2 = Table[pF2[i,j],{j,num},{i,num}];
mS2 =Table[pS2[i,j],{j,num},{i,num}];
SE=mNg2 Map[ListCorrelate[hkern2,#,{2,2},0]&,vars]+mF2  Transpose[Map[ListCorrelate[vkern2,#,{2,2},0]&,Transpose[vars]]]-2 mS2 Transpose[Map[ListCorrelate[vkern1,#,{2,2},0]&,Transpose[Map[ListCorrelate[hkern1,#,{2,2},0]&,vars]]]];
eqns=Thread[Map[Flatten,SE==Table[-2 pQ1[i,j],{j,num},{i,num}]]];
{vec,mat}=CoefficientArrays[eqns,Flatten[vars]];
sol2=Partition[LinearSolve[mat,-vec],num]; //Timing
]

MakePsif:= Module[{},ZeroFrame=Prepend[Append[#1,0],0]&;
sol=Transpose[ZeroFrame/@Transpose[ZeroFrame/@sol2]];
Ψf=ListInterpolation[Transpose[sol],{{-2000.,2000.},{0,1350.}}]]

g8f:=ContourPlot[Ψf[x,z],{x,-1000.`,1000.`},{z,0,1350`},ColorFunction->Hue,DisplayFunction->Identity];
g9f:=VectorPlot[Evaluate[{1000 (∂_({z}))Ψf[x,z],-1000 (∂_({x}))Ψf[x,z]}],{x,-1000.`,1000.`},{z,0.`,1350`},AspectRatio->1,DisplayFunction->Identity];
g11:=Show[GraphicsRow[{Show[{g8f,g9f},FrameLabel->{"x(km)","z(dm)"}],g3}]];
g12:= Show[GraphicsRow[{Show[{g8,g9},FrameLabel->{"x(km)","z(dm)"},PlotLabel->"Quasi-geostrophic"],Show[{g8f,g9f},FrameLabel->{"x(km)","z(dm)"},PlotLabel->"Semi-geostrophic"]}],PlotLabel->"Ageostrophic wind and streamfunction"]

SolveSE;
MakePsif; g12

Cas de la confluence

param=param1;
g4
g5
g7
SolveQGSE;MakePsi;
SolveSE;MakePsif;
g12

Cas de l’advection froide

param=param2;
g4
g5
g7
SolveQGSE;MakePsi;
SolveSE;MakePsif;
g12

Cas de l’advection chaude

param=param3;
g4
g5
g7
SolveQGSE;MakePsi;
SolveSE;MakePsif;
g12

Solution using NDSolve with bad boundary conditions

sol=Flatten[NDSolve[{pNg2[x,z] D[Ψ[z,x],x,x]+(f^2)D[Ψ[z,x],z,z]==- 2 pQ1[x,z] /. param,Ψ[13.5,x]==0 , Derivative[1,0][Ψ][13.5,x]==0 ,Ψ[z,-2000]==Ψ[z,2000] },Ψ,{z,0,13.5},{x,-2000, 2000}]]

ContourPlot[Ψ[z,x] /. sol,{x,-2000,2000},{z,0,13.5},ColorFunction->Hue]

NDSolve[{D[u[x,y],x,x]+D[u[x,y],y,y]==Sin[x]Cos[y],u[0,y]==0, u[2Pi,y]==0,u[x,0]==u[x,2 Pi]},u,{x,0,2 Pi},{y,0, 2 Pi}]

DRAFT AREA

vars=Table[Ψ[x,z],{z,4},{x,4}];

MatrixForm[Map[ListCorrelate[hkern,#,{2,2},0]&,vars]]

MatrixForm[vars]

MatrixForm[Partition[Flatten[vars],4]]

aa=Table[((((i/5)))^2)+0.3 ((((j/5)))^2),{j,50},{i,50}];
ListContourPlot[aa,ColorFunction->Hue]