/[lmdze]/trunk/libf/phylmd/aaam_bud.f90
ViewVC logotype

Contents of /trunk/libf/phylmd/aaam_bud.f90

Parent Directory Parent Directory | Revision Log Revision Log


Revision 61 - (show annotations)
Fri Apr 20 14:58:43 2012 UTC (12 years ago) by guez
File size: 9850 byte(s)
No more included file in LMDZE, not even "netcdf.inc".

Created a variable containing the list of common source files in
GNUmakefile. So we now also see clearly files that are specific to
each program.

Split module "histcom". Assembled resulting files in directory
"Histcom".

Removed aliasing in calls to "laplacien".

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