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

Annotation of /trunk/phylmd/aaam_bud.f

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.21