/[lmdze]/trunk/phylmd/aaam_bud.f
ViewVC logotype

Contents of /trunk/phylmd/aaam_bud.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 82 - (show annotations)
Wed Mar 5 14:57:53 2014 UTC (10 years, 2 months ago) by guez
File size: 9850 byte(s)
Changed all ".f90" suffixes to ".f".
1 module aaam_bud_m
2
3 implicit none
4
5 contains
6
7 subroutine aaam_bud(rea, rg, ome, plat, plon, phis, dragu, liftu, phyu, &
8 dragv, liftv, phyv, p, u, v, aam, torsfc)
9
10 ! Author: F. Lott (LMD/CNRS). Date: 2003/10/20. Object: Compute
11 ! different terms of the axial AAAM budget and mountain torque.
12 ! Only valid for regular rectangular grids. Should be called after
13 ! "lift_noro".
14
15 USE dimens_m, ONLY : iim, jjm
16 use nr_util, only: assert_eq, assert, pi
17
18 real, intent(in):: rea ! Earth radius
19 real, intent(in):: rg ! gravity constant
20 real, intent(in):: ome ! Earth rotation rate
21
22 REAL, intent(in):: plat(:), plon(:)
23 ! (nlon) latitude and longitude in degrees
24
25 real, intent(in):: phis(:) ! (nlon) Geopotential at the ground
26 REAL, intent(in):: dragu(:) ! (nlon) orodrag stress (zonal)
27 REAL, intent(in):: liftu(:) ! (nlon) orolift stress (zonal)
28 REAL, intent(in):: phyu(:) ! (nlon) Stress total de la physique (zonal)
29 REAL, intent(in):: dragv(:) ! (nlon) orodrag stress (Meridional)
30 REAL, intent(in):: liftv(:) ! (nlon) orolift stress (Meridional)
31 REAL, intent(in):: phyv(:) ! (nlon) Stress total de la physique (Meridional)
32
33 REAL, intent(in):: p(:, :)
34 ! (nlon, nlev + 1) pressure (Pa) at model half levels
35
36 real, intent(in):: u(:, :), v(:, :) ! (nlon, nlev) horizontal wind (m/s)
37 REAL, intent(out):: aam ! axial component of wind AAM
38 REAL, intent(out):: torsfc ! axial component of total surface torque
39
40 ! Local Variables:
41
42 INTEGER nlev ! number of vertical levels
43 INTEGER i, j, k, l
44 REAL hadley, hadday
45 REAL dlat, dlon ! latitude and longitude increments (radians)
46
47 REAL raam(3) ! wind AAM (components 1 & 2: equatorial; component 3: axial)
48 REAL oaam(3) ! mass AAM (components 1 & 2: equatorial; component 3: axial)
49 REAL tmou(3) ! resolved mountain torque (3 components)
50 REAL tsso(3) ! parameterised moutain drag torque (3 components)
51 REAL tbls(3) ! parameterised boundary layer torque (3 components)
52 integer iax
53
54 REAL ZS(801, 401) ! topographic height
55 REAL PS(801, 401) ! surface pressure
56 REAL UB(801, 401), VB(801, 401) ! barotropic wind, zonal and meridional
57 REAL SSOU(801, 401), SSOV(801, 401)
58 REAL BLSU(801, 401), BLSV(801, 401)
59 REAL ZLON(801), ZLAT(401) ! longitude and latitude in radians
60
61 !-------------------------------------------------------------------
62
63 call assert(size(plat) == (/size(plon), size(phis), size(dragu), &
64 size(liftu), size(phyu), size(dragv), size(liftv), size(phyv), &
65 size(p, 1), size(u, 1), size(v, 1)/), "aaam_bud nlon")
66 nlev = assert_eq(size(p, 2) - 1, size(u, 2), size(v, 2), "aaam_bud nlev")
67
68 if (iim + 1 > 801 .or. jjm + 1 > 401) then
69 print *, ' Problème de dimension dans aaam_bud'
70 stop 1
71 endif
72
73 hadley = 1e18
74 hadday = 1e18 * 24. * 3600.
75 dlat = pi / jjm
76 dlon = 2 * pi / real(iim)
77
78 oaam = 0.
79 raam = 0.
80 tmou = 0.
81 tsso = 0.
82 tbls = 0.
83
84 ! Mountain height, pressure and barotropic wind:
85
86 ! North pole values (j = 1):
87
88 ub(1, 1) = 0.
89 vb(1, 1) = 0.
90 do k = 1, nlev
91 ub(1, 1) = ub(1, 1) + u(1, k) * (p(1, k) - p(1, k + 1)) / rg
92 vb(1, 1) = vb(1, 1) + v(1, k) * (p(1, k) - p(1, k + 1)) / rg
93 enddo
94
95 zlat(1) = plat(1) * pi / 180.
96
97 do i = 1, iim + 1
98 zs(i, 1) = phis(1) / rg
99 ps(i, 1) = p(1, 1)
100 ub(i, 1) = ub(1, 1)
101 vb(i, 1) = vb(1, 1)
102 ssou(i, 1) = dragu(1) + liftu(1)
103 ssov(i, 1) = dragv(1) + liftv(1)
104 blsu(i, 1) = phyu(1) - dragu(1) - liftu(1)
105 blsv(i, 1) = phyv(1) - dragv(1) - liftv(1)
106 enddo
107
108 l = 1
109 do j = 2, jjm
110 ! Values at Greenwich (Periodicity)
111
112 zs(iim + 1, j) = phis(l + 1) / rg
113 ps(iim + 1, j) = p(l + 1, 1)
114 ssou(iim + 1, j) = dragu(l + 1) + liftu(l + 1)
115 ssov(iim + 1, j) = dragv(l + 1) + liftv(l + 1)
116 blsu(iim + 1, j) = phyu(l + 1) - dragu(l + 1) - liftu(l + 1)
117 blsv(iim + 1, j) = phyv(l + 1) - dragv(l + 1) - liftv(l + 1)
118 zlon(iim + 1) = - plon(l + 1) * pi / 180.
119 zlat(j) = plat(l + 1) * pi / 180.
120
121 ub(iim + 1, j) = 0.
122 vb(iim + 1, j) = 0.
123 do k = 1, nlev
124 ub(iim + 1, j) = ub(iim + 1, j) &
125 + u(l + 1, k) * (p(l + 1, k) - p(l + 1, k + 1)) / rg
126 vb(iim + 1, j) = vb(iim + 1, j) &
127 + v(l + 1, k) * (p(l + 1, k) - p(l + 1, k + 1)) / rg
128 enddo
129
130 do i = 1, iim
131 l = l + 1
132 zs(i, j) = phis(l) / rg
133 ps(i, j) = p(l, 1)
134 ssou(i, j) = dragu(l) + liftu(l)
135 ssov(i, j) = dragv(l) + liftv(l)
136 blsu(i, j) = phyu(l) - dragu(l) - liftu(l)
137 blsv(i, j) = phyv(l) - dragv(l) - liftv(l)
138 zlon(i) = plon(l) * pi / 180.
139
140 ub(i, j) = 0.
141 vb(i, j) = 0.
142 do k = 1, nlev
143 ub(i, j) = ub(i, j) + u(l, k) * (p(l, k) - p(l, k + 1)) / rg
144 vb(i, j) = vb(i, j) + v(l, k) * (p(l, k) - p(l, k + 1)) / rg
145 enddo
146 enddo
147 enddo
148
149 ! South Pole
150
151 l = l + 1
152 ub(1, jjm + 1) = 0.
153 vb(1, jjm + 1) = 0.
154 do k = 1, nlev
155 ub(1, jjm + 1) = ub(1, jjm + 1) + u(l, k) * (p(l, k) - p(l, k + 1)) / rg
156 vb(1, jjm + 1) = vb(1, jjm + 1) + v(l, k) * (p(l, k) - p(l, k + 1)) / rg
157 enddo
158 zlat(jjm + 1) = plat(l) * pi / 180.
159
160 do i = 1, iim + 1
161 zs(i, jjm + 1) = phis(l) / rg
162 ps(i, jjm + 1) = p(l, 1)
163 ssou(i, jjm + 1) = dragu(l) + liftu(l)
164 ssov(i, jjm + 1) = dragv(l) + liftv(l)
165 blsu(i, jjm + 1) = phyu(l) - dragu(l) - liftu(l)
166 blsv(i, jjm + 1) = phyv(l) - dragv(l) - liftv(l)
167 ub(i, jjm + 1) = ub(1, jjm + 1)
168 vb(i, jjm + 1) = vb(1, jjm + 1)
169 enddo
170
171 ! Moment angulaire
172
173 DO j = 1, jjm
174 DO i = 1, iim
175 raam(1) = raam(1) - rea**3 * dlon * dlat * 0.5 * (cos(zlon(i )) &
176 * sin(zlat(j )) * cos(zlat(j )) * ub(i , j ) + cos(zlon(i )) &
177 * sin(zlat(j + 1)) * cos(zlat(j + 1)) * ub(i , j + 1)) &
178 + rea**3 * dlon * dlat * 0.5 * (sin(zlon(i )) * cos(zlat(j )) &
179 * vb(i , j ) + sin(zlon(i )) * cos(zlat(j + 1)) * vb(i , j + 1))
180
181 oaam(1) = oaam(1) - ome * rea**4 * dlon * dlat / rg * 0.5 &
182 * (cos(zlon(i )) * cos(zlat(j ))**2 * sin(zlat(j )) &
183 * ps(i , j ) + cos(zlon(i )) * cos(zlat(j + 1))**2 &
184 * sin(zlat(j + 1)) * ps(i , j + 1))
185
186 raam(2) = raam(2) - rea**3 * dlon * dlat * 0.5 * (sin(zlon(i )) &
187 * sin(zlat(j )) * cos(zlat(j )) * ub(i , j ) + sin(zlon(i )) &
188 * sin(zlat(j + 1)) * cos(zlat(j + 1)) * ub(i , j + 1)) &
189 - rea**3 * dlon * dlat * 0.5 * (cos(zlon(i )) * cos(zlat(j )) &
190 * vb(i , j ) + cos(zlon(i )) * cos(zlat(j + 1)) * vb(i , j + 1))
191
192 oaam(2) = oaam(2) - ome * rea**4 * dlon * dlat / rg * 0.5 &
193 * (sin(zlon(i )) * cos(zlat(j ))**2 * sin(zlat(j )) &
194 * ps(i , j ) + sin(zlon(i )) * cos(zlat(j + 1))**2 &
195 * sin(zlat(j + 1)) * ps(i , j + 1))
196
197 raam(3) = raam(3) + rea**3 * dlon * dlat * 0.5 * (cos(zlat(j))**2 &
198 * ub(i, j) + cos(zlat(j + 1))**2 * ub(i, j + 1))
199
200 oaam(3) = oaam(3) + ome * rea**4 * dlon * dlat / rg * 0.5 &
201 * (cos(zlat(j))**3 * ps(i, j) + cos(zlat(j + 1))**3 &
202 * ps(i, j + 1))
203 ENDDO
204 ENDDO
205
206 ! Couple des montagnes :
207
208 DO j = 1, jjm
209 DO i = 1, iim
210 tmou(1) = tmou(1) - rea**2 * dlon * 0.5 * sin(zlon(i)) &
211 * (zs(i, j) - zs(i, j + 1)) &
212 * (cos(zlat(j + 1)) * ps(i, j + 1) + cos(zlat(j)) * ps(i, j))
213 tmou(2) = tmou(2) + rea**2 * dlon * 0.5 * cos(zlon(i)) &
214 * (zs(i, j) - zs(i, j + 1)) &
215 * (cos(zlat(j + 1)) * ps(i, j + 1) + cos(zlat(j)) * ps(i, j))
216 ENDDO
217 ENDDO
218
219 DO j = 2, jjm
220 DO i = 1, iim
221 tmou(1) = tmou(1) + rea**2 * dlat * 0.5 * sin(zlat(j)) &
222 * (zs(i + 1, j) - zs(i, j)) &
223 * (cos(zlon(i + 1)) * ps(i + 1, j) + cos(zlon(i)) * ps(i, j))
224 tmou(2) = tmou(2) + rea**2 * dlat * 0.5 * sin(zlat(j)) &
225 * (zs(i + 1, j) - zs(i, j)) &
226 * (sin(zlon(i + 1)) * ps(i + 1, j) + sin(zlon(i)) * ps(i, j))
227 tmou(3) = tmou(3) - rea**2 * dlat * 0.5* cos(zlat(j)) &
228 * (zs(i + 1, j) - zs(i, j)) * (ps(i + 1, j) + ps(i, j))
229 ENDDO
230 ENDDO
231
232 ! Couples des differentes friction au sol :
233
234 DO j = 2, jjm
235 DO i = 1, iim
236 tsso(1) = tsso(1) - rea**3 * cos(zlat(j)) * dlon * dlat* &
237 ssou(i, j) * sin(zlat(j)) * cos(zlon(i)) &
238 + rea**3 * cos(zlat(j)) * dlon * dlat* &
239 ssov(i, j) * sin(zlon(i))
240
241 tsso(2) = tsso(2) - rea**3 * cos(zlat(j)) * dlon * dlat* &
242 ssou(i, j) * sin(zlat(j)) * sin(zlon(i)) &
243 - rea**3 * cos(zlat(j)) * dlon * dlat* &
244 ssov(i, j) * cos(zlon(i))
245
246 tsso(3) = tsso(3) + rea**3 * cos(zlat(j)) * dlon * dlat* &
247 ssou(i, j) * cos(zlat(j))
248
249 tbls(1) = tbls(1) - rea**3 * cos(zlat(j)) * dlon * dlat* &
250 blsu(i, j) * sin(zlat(j)) * cos(zlon(i)) &
251 + rea**3 * cos(zlat(j)) * dlon * dlat* &
252 blsv(i, j) * sin(zlon(i))
253
254 tbls(2) = tbls(2) - rea**3 * cos(zlat(j)) * dlon * dlat* &
255 blsu(i, j) * sin(zlat(j)) * sin(zlon(i)) &
256 - rea**3 * cos(zlat(j)) * dlon * dlat* &
257 blsv(i, j) * cos(zlon(i))
258
259 tbls(3) = tbls(3) + rea**3 * cos(zlat(j)) * dlon * dlat* &
260 blsu(i, j) * cos(zlat(j))
261 ENDDO
262 ENDDO
263
264 aam = raam(3)
265 torsfc = tmou(3) + tsso(3) + tbls(3)
266
267 END subroutine aaam_bud
268
269 end module aaam_bud_m

  ViewVC Help
Powered by ViewVC 1.1.21