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

Annotation of /trunk/Sources/phylmd/aaam_bud.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 76 - (hide annotations)
Fri Nov 15 18:45:49 2013 UTC (10 years, 6 months ago) by guez
Original Path: trunk/phylmd/aaam_bud.f90
File size: 9850 byte(s)
Moved everything out of libf.
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     integer iax
53 guez 3
54 guez 56 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 guez 3
61 guez 56 !-------------------------------------------------------------------
62 guez 3
63 guez 56 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 guez 3
68 guez 56 if (iim + 1 > 801 .or. jjm + 1 > 401) then
69     print *, ' Problème de dimension dans aaam_bud'
70     stop 1
71     endif
72 guez 3
73 guez 56 hadley = 1e18
74     hadday = 1e18 * 24. * 3600.
75     dlat = pi / jjm
76     dlon = 2 * pi / real(iim)
77 guez 3
78 guez 56 oaam = 0.
79     raam = 0.
80     tmou = 0.
81     tsso = 0.
82     tbls = 0.
83 guez 3
84 guez 56 ! Mountain height, pressure and barotropic wind:
85 guez 3
86 guez 56 ! North pole values (j = 1):
87 guez 3
88 guez 56 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 guez 3
95 guez 56 zlat(1) = plat(1) * pi / 180.
96 guez 3
97 guez 56 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 guez 3
108 guez 56 l = 1
109     do j = 2, jjm
110     ! Values at Greenwich (Periodicity)
111 guez 3
112 guez 56 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 guez 3
121 guez 56 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 guez 3
130 guez 56 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 guez 3
140 guez 56 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 guez 3
149 guez 56 ! South Pole
150 guez 3
151 guez 56 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 guez 3
160 guez 56 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 guez 3
171 guez 56 ! Moment angulaire
172 guez 3
173 guez 56 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 guez 3
181 guez 56 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 guez 3
186 guez 56 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 guez 3
192 guez 56 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 guez 3
197 guez 56 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 guez 3
200 guez 56 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 guez 3
206 guez 56 ! Couple des montagnes :
207 guez 3
208 guez 56 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 guez 3
219 guez 56 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 guez 3
232 guez 56 ! 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