-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathBtransfo.f90
165 lines (160 loc) · 4.39 KB
/
Btransfo.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
Subroutine q2psi(psi,q)
!c
!c Transformation Q-->PSI
!c
include 'sphectra.h'
include 'paramod.h'
!c
real psi(nvaria2,nlevels),q(nvaria2,nlevels)
real trv(nvaria2,nlevels)
real trv1(nlevels),trv2(nlevels)
!c
!c Translation de q pour obtenir le TP relatif
!c C est Q - f - fh/Ho
!c
call scopy(nvaria2*nlevels,q,1,trv,1)
call saxpy(nvaria2,-1.,ftopo,1,trv(1,3),1)
do nl=1,nlevels
trv(1,nl) = trv(1,nl) - f1
enddo
!c
!c Partie lineaire de la transformation
!c
do nv=1,nvaria2
do nl=1,nlevels
trv1(nl) = trv(nv,nl)
enddo
call murrv(nlevels,nlevels,cq2psi(1,1,nv),nlevels,nlevels,trv1,1,nlevels,trv2)
do nl=1,nlevels
psi(nv,nl) = trv2(nl)
enddo
enddo
!c
return
end
!c
subroutine psi2q(q,psi)
!c
!c Programme passant de psi, la fonction de courant (geost.), a q,
!c la vorticite potentielle.
!c
include 'sphectra.h'
include 'paramod.h'
!c
real psi(nvaria2,nlevels),q(nvaria2,nlevels)
!c
call sset(nvaria2*nlevels,0.,q,1)
!c f1=omega2/sqrt(3.)
!c namax=2*mvadeb(ntrunc-1)
!c
!c
do nv2=1,nvaria2
oplap = oprlap((nv2-1)/2+1,2)
q(nv2,1)=oplap*psi(nv2,1)-(psi(nv2,1)-psi(nv2,2))*ray2n2
q(nv2,2)=oplap*psi(nv2,2) &
-(psi(nv2,2)-psi(nv2,3))*ray4n2 &
+(psi(nv2,1)-psi(nv2,2))*ray2n2
q(nv2,3)=oplap*psi(nv2,3) &
+(psi(nv2,2)-psi(nv2,3))*ray4n2 &
+ftopo(nv2)
enddo
!c
!c Addition de La vorticite planetaire
!c
do nl=1,nlevels
q(1,nl)=q(1,nl)+f1
enddo
!c
return
end
!c
subroutine dq2dpsi(psi,q)
!c
!c Transformation Q-->PSI. Ne sont retenus que les termes lineaires
!c Les termes constants sont retires.
!c
include 'sphectra.h'
include 'paramod.h'
!c
real psi(nvaria2,nlevels),q(nvaria2,nlevels)
real trv(nvaria2,nlevels)
real trv1(nlevels),trv2(nlevels)
!c
call scopy(nvaria2*nlevels,q,1,trv,1)
!c
!c Partie lineaire de la transformation
!c
do nv=1,nvaria2
do nl=1,nlevels
trv1(nl) = trv(nv,nl)
enddo
call murrv(nlevels,nlevels,cq2psi(1,1,nv),nlevels,nlevels,trv1,1,nlevels,trv2)
do nl=1,nlevels
psi(nv,nl) = trv2(nl)
enddo
enddo
!c
return
end
!c
subroutine dpsi2dq(q,psi)
!c
!c Programme passant de psi, la fonction de courant (geost.), a q,
!c la vorticite potentielle. N est retenu ici que la partie lineaire
!c de la transformation
!c
include 'sphectra.h'
include 'paramod.h'
!c
real psi(nvaria2,nlevels),q(nvaria2,nlevels)
!c
call sset(nvaria2*nlevels,0.,q,1)
!c f1=omega2/sqrt(3.)
!c namax=2*mvadeb(ntrunc-1)
!c
!c
do nv2=1,nvaria2
oplap = oprlap((nv2-1)/2+1,2)
q(nv2,1)=oplap*psi(nv2,1)-(psi(nv2,1)-psi(nv2,2))*ray2n2
q(nv2,2)=oplap*psi(nv2,2) &
-(psi(nv2,2)-psi(nv2,3))*ray4n2 &
+(psi(nv2,1)-psi(nv2,2))*ray2n2
q(nv2,3)=oplap*psi(nv2,3) &
+(psi(nv2,2)-psi(nv2,3))*ray4n2
enddo
!c
return
end
!c
subroutine zetap2q(q,zetap)
!c On a q=lap(psi)+R psi + b = zetap +b. Le vecteur b est le vecteur constant
!c contennant la vorticite planetaire et la topographie.
include 'sphectra.h'
include 'paramod.h'
real zetap(nvaria2,nlevels),q(nvaria2,nlevels)
do nl=1,nlevels
call scopy(nvaria2,zetap(1,nl),1,q(1,nl),1)
q(1,nl) = q(1,nl) + f1
if(nl.eq.nlevels) then
call saxpy(nvaria2,1.,ftopo,1,q(1,nl),1)
endif
enddo
return
end
!c
!c
subroutine q2zetap(zetap,q)
!c On a q=lap(psi)+R psi + b = zetap +b. Le vecteur b est le vecteur constant
!c contennant la vorticite planetaire et la topographie.
include 'sphectra.h'
include 'paramod.h'
real zetap(nvaria2,nlevels),q(nvaria2,nlevels)
do nl=1,nlevels
call scopy(nvaria2,q(1,nl),1,zetap(1,nl),1)
zetap(1,nl) = zetap(1,nl) - f1
if(nl.eq.nlevels) then
call saxpy(nvaria2,-1.,ftopo,1,zetap(1,nl),1)
endif
enddo
return
end