Bennett A.F. Inverse Modeling of the Ocean and Atmosphere
Подождите немного. Документ загружается.
A.3 Discrete formulation 207
(λ
v
)
k+1
i, j
≡ W
v
f
ˆv
k+1
i, j
− ˆv
k
i, j
t
+ f
ˆ
u
k
i+1, j
+
ˆ
u
k
i, j
+
ˆ
u
k
i+1, j −1
+
ˆ
u
k
i, j −1
4
+ g
ˆ
q
k+1
i, j
−
ˆ
q
k+1
i, j −1
y
+r
v
ˆv
k
i, j
− (F
v
)
k
i, j
,
(λ
q
)
k+1
i, j
≡ W
q
f
ˆ
q
k+1
i, j
−
ˆ
q
k
i, j
t
+ H
ˆ
u
k
i+1, j
−
ˆ
u
k
i, j
x
+
ˆv
k
i, j +1
− ˆv
k
i, j
y
+r
q
ˆ
q
k
i, j
.
Euler–Lagrange equations
−
(λ
u
)
k+1
i, j
− (λ
u
)
k
i, j
t
+ f
(λ
v
)
k+1
i, j +1
+ (λ
v
)
k+1
i, j
+ (λ
v
)
k+1
i−1, j +1
+ (λ
v
)
k+1
i−1, j
4
− H
(λ
q
)
k+1
i, j
− (λ
q
)
k+1
i−1, j
x
+r
u
(λ
u
)
k+1
i, j
= 0,
−
(λ
v
)
k+1
i, j
− (λ
v
)
k
i, j
t
− f
(λ
u
)
k+1
i+1, j
+ (λ
u
)
k+1
i, j
+ (λ
u
)
k+1
i+1, j −1
+ (λ
u
)
k+1
i, j −1
4
− H
(λ
q
)
k+1
i, j
− (λ
q
)
k+1
i, j −1
y
+r
v
(λ
v
)
k+1
i, j
= 0,
−
(λ
q
)
k+1
i, j
− (λ
q
)
k
i, j
t
− g
(λ
u
)
k
i+1, j
− (λ
u
)
k
i, j
x
+
(λ
v
)
k
i, j +1
− (λ
v
)
k
i, j
y
+ r
q
(λ
q
)
k+1
i, j
=−w
M
m=1
δ
i,i
m
δ
j, j
m
δ
k,k
m
xyt
ˆ
q
k
i
m
, j
m
− d
m
,
where k = 1, NK − 1.
(λ
u
)
NK
i, j
t
= 0,
(λ
v
)
NK
i, j
t
= 0,
(λ
q
)
NK
i, j
t
− g
(λ
u
)
NK
i+1, j
− (λ
u
)
NK
i, j
x
+
(λ
v
)
NK
i, j +1
− (λ
v
)
NK
i, j
y
=−w
M
m=1
δ
i,i
m
δ
j, j
m
δ
NK,k
m
xyt
ˆ
q
NK
i
m
, j
m
− d
m
.
(λ
v
)
k
i,1
= 0,
(λ
v
)
k
i,NJ
= 0 (computational boundary conditions).
208 Appendix A Computing exercises
(λ
u
)
k
0, j
= (λ
u
)
k
NI, j
,
(λ
u
)
k
NI+1, j
= (λ
u
)
k
1, j
,
(λ
v
)
k
0, j
= (λ
v
)
k
NI, j
,
(λ
v
)
k
NI+1, j
= (λ
v
)
k
1, j
,
(λ
q
)
k
0, j
= (λ
q
)
k
NI, j
,
(λ
q
)
k
NI+1, j
= (λ
q
)
k
1, j
.
ˆ
u
k+1
i, j
−
ˆ
u
k
i, j
t
− f
ˆv
k
i, j +1
+ ˆv
k
i, j
+ ˆv
k
i−1, j +1
+ ˆv
k
i−1, j
4
+ g
ˆ
q
k+1
i, j
−
ˆ
q
k+1
i−1, j
x
+r
u
ˆ
u
k
i, j
= (F
u
)
k
i, j
+
'
W
u
f
(
−1
(λ
u
)
k+1
i, j
,
ˆv
k+1
i, j
− ˆv
k
i, j
t
+ f
ˆ
u
k
i+1, j
+
ˆ
u
k
i, j
+
ˆ
u
k
i+1, j −1
+
ˆ
u
k
i, j −1
4
+ g
ˆ
q
k+1
i, j
−
ˆ
q
k+1
i, j −1
y
+r
v
ˆv
k
i, j
= (F
v
)
k
i, j
+
'
W
v
f
(
−1
(λ
v
)
k+1
i, j
,
ˆ
q
k+1
i, j
−
ˆ
q
k
i, j
t
+ H
ˆ
u
k
i+1, j
−
ˆ
u
k
i, j
x
+
ˆv
k
i, j +1
− ˆv
k
i, j
y
+r
q
ˆ
q
k
i, j
=
'
W
q
f
(
−1
(λ
q
)
k+1
i, j
.
W
u
i
ˆ
u
0
i, j
− (I
u
)
i, j
t
=
(λ
u
)
1
i, j
t
− f
(λ
v
)
1
i, j +1
+ (λ
v
)
1
i, j
+ (λ
v
)
1
i−1, j +1
+ (λ
v
)
1
i−1, j
4
+ H
(λ
q
)
1
i, j
− (λ
q
)
1
i−1, j
x
,
W
v
i
ˆv
0
i, j
− (I
v
)
i, j
t
=
(λ
v
)
1
i, j
t
+ f
(λ
u
)
1
i+1, j
+ (λ
u
)
1
i, j
+ (λ
u
)
1
i+1, j −1
+ (λ
u
)
1
i, j −1
4
+ H
(λ
q
)
1
i, j
− (λ
q
)
1
i, j −1
y
−r
v
(λ
v
)
1
i, j
,
W
q
i
ˆ
q
0
i, j
− (I
q
)
i, j
t
=
(λ
q
)
1
i, j
t
−r
q
(λ
q
)
1
i, j
.
ˆv
k
i,1
=
'
W
v
b
(
−1
y
f
(λ
u
)
k
i,1
+ (λ
u
)
k
i+1,1
4
+ H
(λ
q
)
k+1
i,1
y
,
ˆv
k
i,NJ
=
'
W
v
b
(
−1
y
f
(λ
u
)
k
i,NJ−1
+ (λ
u
)
k
i+1,NJ−1
4
− H
(λ
q
)
k+1
i,NJ−1
y
.
ˆ
u
k
0, j
=
ˆ
u
k
NI, j
,
ˆ
u
k
NI+1, j
=
ˆ
u
k
1, j
,
A.4 Representer calculation 209
ˆv
k
0, j
= ˆv
k
NI, j
,
ˆv
k
NI+1, j
= ˆv
k
1, j
,
ˆ
q
k
0, j
=
ˆ
q
k
NI, j
,
ˆ
q
k
NI+1, j
=
ˆ
q
k
1, j
.
A.4
Representer calculation
A.4.1 Preamble
The first part of the exercise is to calculate some representers, assemble the representer
matrix and verify its algebraic properties. The second part is to find the extremum of a
penalty functional by solving the EL equations.
A.4.2
Representer vector
Calculate the representer vector r(x, y, t).
Source code: rep.f
To compile the source code: f77-O3 rep.f -o rep
A.4.3
Representer matrix
Construct the representer matrix R.
Check: Is R symmetric and positive-definite?
A.4.4
Extremum
Find the extremum of the penalty functional J .
Note: Avoid storing r(x, y, t ) when assembling (3.24).
Instead, substitute for the coupling, integrate the backward EL equations and
then integrate the forward equations (see §3.1: Accelerating the representer
calculation).
Check: Is
ˆ
β
m
=−w
[
ˆ
q(x
m
, y
m
, t
m
) − d
m
]
?
A.4.5
Weights
The following values are suggested:
Dynamical weight (u)
W
u
f
xyt = (0.25|F
u
|)
−2
s
4
m
−2
.
210 Appendix A Computing exercises
Dynamical weight (v)
W
v
f
xyt = (0.25|F
u
|)
−2
s
4
m
−2
.
Dynamical weight (q)
W
q
f
xyt =∞.
Initial weight (u)
W
u
i
xy =∞.
Initial weight (v)
W
v
i
xy =∞.
Initial weight (q)
W
q
i
xy =∞.
Boundary weight (v)
W
v
b
xt =∞.
Data weight (q)
w = [0.1 max(q
F
)]
−2
m
−2
.
A.5
More representer calculations
A.5.1 Pseudocode, preconditioned conjugate gradient solver
Find the local extremum of the penalty functional J iteratively (see §3.1, Accelerating
the representer calculation). Below is the “pseudocode” of a preconditioned conjugate
gradient method (Golub & Van Loan, 1989) for solving
U
−1/2
(R + w
−1
I)U
−1/2
U
1/2
ˆ
β = U
−1/2
h,
where U is a preconditioner.
l = 0;
ˆ
β
0
= 0; e
0
= h; ω
0
=!e
0
!
2
/! h !
2
while ω
l
<ω
sc
solve z
l
= U
−1
e
l
l = l + 1
if l = 1
p
1
= z
0
else
γ
l
= e
T
l−1
z
l−1
/e
T
l−2
z
l−2
p
l
= z
l−1
+ γ
l
p
l−1
endif
α
l
= e
T
l−1
z
l−1
/p
T
l
(R + w
−1
I)p
l
ˆ
β
l
=
ˆ
β
l−1
+ α
l
p
l
A.5 More representer calculations 211
e
l
= e
l−1
− α
l
(R + w
−1
I)p
l
ω
l
=!e
l
!
2
/! h !
2
endwhile
ˆ
β =
ˆ
β
l
Source code: cgm.f (conjugate gradient solver).
A.5.2
Convolutions for covarying errors
Use “non-diagonal” covariances for initial and dynamical errors (see §2.6, Smoothing
norms, covariances and convolutions).
Appendix B
Euler–Lagrange equations for a numerical
weather prediction model
The dynamics are those of the standard σ -coordinate, Primitive-Equation model of
a moist atmosphere on the sphere (Haltiner and Williams, 1980, p. 17). A penalty
functional and the associated Euler–Lagrange equations are given in continuous form;
CMFortran code for finite-difference forms is available at an anonymous ftp site.
Details of the measurement functionals for reprocessed cloud-track wind observations
(see §5.4), and the associated impulses in the adjoint equations, have been suppressed
here. The details may be found in the code.
B.1
Symbols
a
0
," earth’s radius, rotation rate
λ, φ, σ, t longitude, latitude, sigma, time
u,v, ˙σ zonal, meridional, vertical velocity components
f ≡ 2" sin φ Coriolis parameter
#, T, T
v
, p, p
∗
= p/σ geopotential, temperature, virtual temperature, pressure,
surface pressure
q,ρ specific humidity, density
l ≡ ln p
∗
log surface pressure
R
d
, C
pd
gas constant, specific heat at constant pressure
(both for dry air)
R
v
, C
pv
gas constant, specific heat at constant pressure
(both for water vapor)
= R
v
/R
d
,δ= C
pv
/C
pd
212
B.2 Primitive Equations and penalty functional 213
S
u
, S
v
, S
T
, S
q
, S
l
prior estimates of sources of zonal momentum,
meridional momentum, heat, humidity, mass
u≡
1
!
0
udσ, u
≡ u −u vertical average, fluctuation
ρ
u
,ρ
v
,ρ
T
,ρ
q
,ρ
l
,ρ
σ
dynamical residuals, or errors in prior estimates
of sources
W
u
, W
T
, W
q
, W
#
, W
l
, W
σ
weights for residuals (inverses of prior estimates
of covariances of dynamical residuals)
Q
u
, Q
T
, Q
q
, Q
#
, Q
l
, Q
σ
prior estimates of covariances of dynamical
residuals
µ, θ, κ, χ, ξ, ω weighted residuals, or adjoint variables
J = J [u, ˙σ,#,T, q, l] penalty functional, or estimator for residuals
F ≡ ∂ T /∂T
v
= 1 + (
−1
− 1)q moisture factor
D ≡
1
a
0
cos θ
u
λ
+
(v cos θ )
θ
a
0
cos θ
+
∂ ˙σ
∂σ
−
˙σ
σ
divergence
[δu], [δv], [δ#], [δT ], Euler–Lagrange equations for arbitrary variations
[δ ˙σ ], [δl], [δq]ofu,v,#,T, ˙σ,l, q
u
I
, T
I
, q
I
, l
I
prior estimates of initial values for u, T, q, l
V
u
, V
T
, V
q
, V
l
weights for residuals or errors in prior estimates of
initial values (inverses of initial error covariances)
O
u
, O
T
, O
q
, O
l
prior estimates of initial error covariances for
u, T, q, l
#
∗
, V
∗
, O
∗
orography, weight, error covariance
•, ◦ four-dimensional and three-dimensional inner
products
r
n
u
, r
n
v
, r
n
T
, r
n
#
, r
n
q
, r
n
l
, r
n
˙σ
representers for n
th
iterate of Euler–Lagrange
equations
a
n
u
, a
n
v
, a
n
T
, a
n
φ
, a
n
q
, a
n
l
, a
n
˙σ
adjoint representers
B.2
Primitive Equations and penalty functional
u
t
+
u
a
0
cos φ
u
λ
+
v
a
0
u
φ
+ ˙σ u
σ
−
f +
u
a
0
tan φ
v
+
1
a
0
cos φ
(#
λ
+ R
d
T
v
l
λ
) − S
u
≡ ρ
u
(B.2.1)
v
t
+
u
a
0
cos φ
v
λ
+
v
a
0
v
φ
+ ˙σv
σ
+
f +
u
a
0
tan φ
u
+
1
a
0
(#
φ
+ R
d
T
v
l
φ
) − S
v
≡ ρ
v
(B.2.2)
214 Appendix B Numerical weather prediction
T
t
+
u
a
0
cos φ
T
λ
+
v
a
0
T
φ
+ ˙σ T
σ
+
R
d
T
v
C
pd
B
1
a
0
cos φ
u
λ
+
(v cos φ)
φ
a
0
cos φ
+σ
˙σ
σ
σ
− S
T
≡ ρ
T
(B.2.3)
q
t
+
u
a
0
cos φ
q
λ
+
v
a
0
q
φ
+ ˙σ q
σ
− S
q
≡ ρ
q
∂#
∂lnσ
+ R
d
T
v
≡ ρ
#
(B.2.4)
l
t
+
u
a
0
cos φ
l
λ
+
v
a
0
l
φ
+
1
a
0
cos φ
u
λ
+
(vcos φ)
φ
a
0
cos φ
−S
l
≡ρ
l
(B.2.5)
u
a
0
cos φ
l
λ
+
v
a
0
l
φ
+
1
a
0
cos φ
u
λ
+
(v
cos φ)
φ
a
0
cos φ
+ ˙σ
σ
− S
l
≡ ρ
σ
(B.2.6)
T
v
≡ T [1 + (ε
−1
− 1)q], p = ρ R
d
T
v
, B ≡ 1 +(δ − 1)q. (B.2.7)
Weighted residuals
µ ≡ W
u
• ρ
u
,θ≡ W
T
• ρ
T
,κ≡ W
q
• ρ
q
, (B.2.8)
χ = W
#
• ρ
#
,ξ= W
l
• ρ
l
,ω= W
σ
• ρ
σ
. (B.2.9)
Penalty functional
J = J [u, ˙σ,#,T, q, l] = ρ
∗
u
• W
u
• ρ
u
+ ρ
T
• W
T
• ρ
T
+ ρ
q
• W
q
• ρ
q
+ρ
#
• W
#
• ρ
#
+ ρ
l
• W
l
• ρ
l
+ ρ
σ
• W
σ
• ρ
σ
+ (boundary penalties @ σ = 0, 1) +(initial penalties for u, T, q, & l)
+ (data penalties). (B.2.10)
B.3
Euler–Lagrange equations
Note: The symbols [δu], etc., indicate that the following equation is the extremal
condition for the penalty functional J , with respect to variations of δu, etc.
[δu]
−µ
t
+
u
λ
µ
a
0
cos φ
−
(uµ)
λ
a
0
cos φ
−
(vµ cos φ)
φ
a
0
cos φ
− (˙σµ)
σ
−
tan φ
a
0
vµ +
v
λ
ν
a
0
cos φ
+
f +
u tan φ
a
0
ν +
tan φ
a
0
uν
+
T
λ
θ
a
0
cos φ
−
R
d
C
pd
T
v
θ
B
λ
1
a
0
cos φ
+
q
λ
κ
a
0
cos φ
+
l
λ
ξ
a
0
cos φ
−
ξ
λ
a
0
cos φ
+
l
λ
(ω −ω)
a
0
cos φ
−
(ω
λ
−ω
λ
)
a
0
cos φ
− (impulses) = 0 (B.3.1)
B.3 Euler–Lagrange equations 215
[δv]
−ν
t
+
u
φ
µ
a
0
−
f +
u tan φ
a
0
µ −
(uν)
λ
a
0
cos φ
+
v
φ
ν
a
0
−
(vν cos φ)
φ
a
0
cos φ
− (˙σν)
σ
+
T
φ
θ
a
0
−
R
d
C
pd
T
v
θ
B
φ
1
a
0
+
q
φ
κ
a
0
+
l
φ
ξ
a
0
−
ξ
φ
a
0
+
l
φ
(ω −ω)
a
0
−
(ω
φ
−ω
φ
)
a
0
− (impulses) = 0 (B.3.2)
[δ#]
−
µ
λ
a
0
cos φ
−
(ν cos φ)
φ
a
0
cos φ
− (σχ)
σ
− (impulses) = 0 (B.3.3)
[δT ]
−θ
t
+
R
d
Fl
λ
µ
a
0
cos φ
+
R
d
Fl
φ
ν
a
0
−
(uθ )
λ
a
0
cos φ
−
(vθ cos φ)
φ
a
0
cos φ
− (˙σθ)
σ
+
R
d
FDθ
C
pd
B
+ R
d
Fχ − (impulses) = 0 (B.3.4)
[δ ˙σ ]
u
σ
µ + v
σ
ν + T
σ
θ −
R
d
C
pd
T
v
θ
B
σ
−
R
d
C
pd
T
v
θ
Bσ
+ q
σ
κ −ω
σ
− (impulses) = 0
(B.3.5)
[δl]
−ξ
t
−
R
d
(T
v
µ)
λ
a
0
cos φ
−
R
d
(T
v
ν cos φ)
φ
a
0
cos φ
−
(uξ )
λ
a
0
cos φ
−
(vξ cos φ)
φ
a
0
cos φ
−
u
ω
λ
a
0
cos φ
−
v
ω cos φ
φ
a
0
cos φ
− (impulses) = 0 (B.3.6)
[δq]
−κ
t
+
R
d
TF
q
l
λ
µ
a
0
cos φ
+
R
d
TF
q
l
φ
ν
a
0
+
R
d
T
C
pd
F
q
Dθ
B
−
R
d
C
pd
Dθ(δ − 1)T
v
B
2
−
(uκ)
λ
a
0
cos φ
−
(vκ cos φ)
φ
a
0
cos φ
+ R
d
F
q
χ −(˙σκ)
σ
−(impulses) = 0 (B.3.7)
Note: (B.3.5) is a first-order, ordinary differential equation for ω as a function of σ .
Solutions are indeterminate, since there is no boundary condition for ω at
σ = 0, nor at σ = 1. However, the other Euler–Lagrange equations only
involve ω
or u
ω, where ω
≡ ω −ω, etc., thus the indeterminacy has no
effect upon them. The residual ρ
σ
in (B.2.6) must satisfy ρ
σ
=0; thus its
covariance Q
σ
must have vanishing integrals with respect to both vertical
216 Appendix B Numerical weather prediction
arguments, and any hypothesis for Q
σ
must conform to this requirement. Hence
the optimal estimate ρ
σ
≡ Q
σ
• ω is also unaffected by the indeterminacy, and
the vertical integral of the estimate vanishes.
Initial conditions
@ t = 0: u
∼
=
u
I
, T
∼
=
T
I
, q
∼
=
q
I
, l
∼
=
l
I
. (B.3.8)
Contribution to penalty functional
J
I
= (u − u
I
)
∗
◦ V
u
◦ (u − u
I
) + (T − T
I
) ◦ V
T
◦ (T − T
I
)
+ (q − q
I
) ◦ V
q
◦ (q − q
I
) + (l −l
I
) ◦ V
l
◦ (l − l
I
). (B.3.9)
Hence
@ t = T : µ = 0 ,θ= 0 ,κ= 0 ,ξ= 0 (B.3.10)
@t = 0: −µ + V
u
◦ (u − u
I
) = 0 , −θ + V
T
◦ (T − T
I
) = 0,
−κ + V
q
◦ (q − q
I
) = 0, −ξ + V
l
◦ (l − l
I
) = 0, (B.3.11)
i.e., u = u
I
+ O
u
◦ µ, T = T
I
+ O
T
◦ θ, q = q
I
+ O
q
◦ κ, l = l
I
+ O
l
◦ ξ,
where O
u
(12)
◦ V
u
(23)
= δ(x
1
− x
3
)I, etc. (B.3.12)
Boundary conditions
@ σ = 0, 1: ˙σ = 0 (B.3.13)
@ σ = 1: #
∼
=
#
∗
. (B.3.14)
Contribution to penalty functional
J
∗
= (# − #
∗
) ◦ V
∗
◦ (# − #
∗
). (B.3.15)
Hence
@ σ = 0: χ = 0, (B.3.16)
@ σ = 1: # = #
∗
− O
∗
◦ χ. (B.3.17)
B.4
Linearized Primitive Equations
Note: The labels (LPE1) etc. refer to lines of the CMFortran code for the
finite-difference equations; the code is available at an anonymous ftp site
(ftp.oce.orst.edu, /dist/chua/IOM/IOSU).
u
n
t
+
u
n−1
u
n
λ
a
0
cos φ
+
v
n−1
u
n
φ
a
0
+ ˙σ
n−1
u
n
σ
−
f +
u
n−1
tan φ
a
0
v
n
+
1
a
0
cos φ
#
n
λ
+ R
d
T
n−1
v
l
n
λ
− S
n
u
≡ ρ
n
u
. (LPE1) (B.4.1)