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

Contents of /trunk/phylmd/aaam_bud.f

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.21