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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 61 - (hide annotations)
Fri Apr 20 14:58:43 2012 UTC (12 years, 1 month 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 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