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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 47 - (show annotations)
Fri Jul 1 15:00:48 2011 UTC (12 years, 9 months ago) by guez
File size: 10737 byte(s)
Split "thermcell.f" and "cv3_routines.f".
Removed copies of files that are now in "L_util".
Moved "mva9" and "diagetpq" to their own files.
Unified variable names across procedures.

1 subroutine aaam_bud (iam,nlon,nlev,rsec,
2 i rea,rg,ome,
3 i plat,plon,phis,
4 i dragu,liftu,phyu,
5 i dragv,liftv,phyv,
6 i p, u, v,
7 o aam, torsfc)
8 c
9 use dimens_m
10 use dimphy
11 implicit none
12 c======================================================================
13 c Auteur(s): F.Lott (LMD/CNRS) date: 20031020
14 c Object: Compute different terms of the axial AAAM Budget.
15 C No outputs, every AAM quantities are written on the IAM
16 C File.
17 c
18 c Modif : I.Musat (LMD/CNRS) date : 20041020
19 c Outputs : axial components of wind AAM "aam" and total surface torque
20 C "torsfc",
21 c but no write in the iam file.
22 c
23 C WARNING: Only valid for regular rectangular grids.
24 C REMARK: CALL DANS PHYSIQ AFTER lift_noro:
25 C CALL aaam_bud (27,klon,klev, gmtime,
26 C C ra,rg,romega,
27 C C rlat,rlon,pphis,
28 C C zustrdr,zustrli,zustrph,
29 C C zvstrdr,zvstrli,zvstrph,
30 C C paprs,u,v)
31 C
32 C======================================================================
33 c Explicit Arguments:
34 c ==================
35 c iam-----input-I-File number where AAMs and torques are written
36 c It is a formatted file that has been opened
37 c in physiq.F
38 c nlon----input-I-Total number of horizontal points that get into physics
39 c nlev----input-I-Number of vertical levels
40 c rsec----input-R-Seconde de la journee
41 c rea-----input-R-Earth radius
42 c rg------input-R-gravity constant
43 c ome-----input-R-Earth rotation rate
44 c plat ---input-R-Latitude en degres
45 c plon ---input-R-Longitude en degres
46 c phis ---input-R-Geopotential at the ground
47 c dragu---input-R-orodrag stress (zonal)
48 c liftu---input-R-orolift stress (zonal)
49 c phyu----input-R-Stress total de la physique (zonal)
50 c dragv---input-R-orodrag stress (Meridional)
51 c liftv---input-R-orolift stress (Meridional)
52 c phyv----input-R-Stress total de la physique (Meridional)
53 c p-------input-R-Pressure (Pa) at model half levels
54 c u-------input-R-Horizontal wind (m/s)
55 c v-------input-R-Meridional wind (m/s)
56 c aam-----output-R-Axial Wind AAM (=raam(3))
57 c torsfc--output-R-Total surface torque (=tmou(3)+tsso(3)+tbls(3))
58 c
59 c Implicit Arguments:
60 c ===================
61 c
62 c iim--common-I: Number of longitude intervals
63 c jjm--common-I: Number of latitude intervals
64 c klon-common-I: Number of points seen by the physics
65 c iim*(jjm-1)+2 for instance
66 c klev-common-I: Number of vertical layers
67 c======================================================================
68 c Local Variables:
69 c ================
70 c dlat-----R: Latitude increment (Radians)
71 c dlon-----R: Longitude increment (Radians)
72 c raam ---R: Wind AAM (3 Components, 1 & 2 Equatoriales; 3 Axiale)
73 c oaam ---R: Mass AAM (3 Components, 1 & 2 Equatoriales; 3 Axiale)
74 c tmou-----R: Resolved Mountain torque (3 components)
75 c tsso-----R: Parameterised Moutain drag torque (3 components)
76 c tbls-----R: Parameterised Boundary layer torque (3 components)
77 c
78 c LOCAL ARRAY:
79 c ===========
80 c zs ---R: Topographic height
81 c ps ---R: Surface Pressure
82 c ub ---R: Barotropic wind zonal
83 c vb ---R: Barotropic wind meridional
84 c zlat ---R: Latitude in radians
85 c zlon ---R: Longitude in radians
86 c======================================================================
87
88 c
89 c ARGUMENTS
90 c
91 INTEGER iam,nlon,nlev
92 real, intent(in):: rsec
93 real rea
94 real, intent(in):: rg
95 real ome
96 REAL, intent(in):: plat(nlon),plon(nlon)
97 real phis(nlon)
98 REAL dragu(nlon),liftu(nlon),phyu(nlon)
99 REAL dragv(nlon),liftv(nlon),phyv(nlon)
100 REAL, intent(in):: p(nlon,nlev+1)
101 real, intent(in):: u(nlon,nlev), v(nlon,nlev)
102 c
103 c Variables locales:
104 c
105 INTEGER i,j,k,l
106 REAL xpi,hadley,hadday
107 REAL dlat,dlon
108 REAL raam(3),oaam(3),tmou(3),tsso(3),tbls(3)
109 integer iax
110 cIM ajout aam, torsfc
111 c aam = composante axiale du Wind AAM raam
112 c torsfc = composante axiale de (tmou+tsso+tbls)
113 REAL aam, torsfc
114
115 REAL ZS(801,401),PS(801,401)
116 REAL UB(801,401),VB(801,401)
117 REAL SSOU(801,401),SSOV(801,401)
118 REAL BLSU(801,401),BLSV(801,401)
119 REAL ZLON(801),ZLAT(401)
120 C
121 C PUT AAM QUANTITIES AT ZERO:
122 C
123 if(iim+1.gt.801.or.jjm+1.gt.401)then
124 print *,' Pb de dimension dans aaam_bud'
125 stop
126 endif
127
128 xpi=acos(-1.)
129 hadley=1.e18
130 hadday=1.e18*24.*3600.
131 dlat=xpi/float(jjm)
132 dlon=2.*xpi/float(iim)
133
134 do iax=1,3
135 oaam(iax)=0.
136 raam(iax)=0.
137 tmou(iax)=0.
138 tsso(iax)=0.
139 tbls(iax)=0.
140 enddo
141
142 C MOUNTAIN HEIGHT, PRESSURE AND BAROTROPIC WIND:
143
144 C North pole values (j=1):
145
146 l=1
147
148 ub(1,1)=0.
149 vb(1,1)=0.
150 do k=1,nlev
151 ub(1,1)=ub(1,1)+u(l,k)*(p(l,k)-p(l,k+1))/rg
152 vb(1,1)=vb(1,1)+v(l,k)*(p(l,k)-p(l,k+1))/rg
153 enddo
154
155 zlat(1)=plat(l)*xpi/180.
156
157 do i=1,iim+1
158
159 zs(i,1)=phis(l)/rg
160 ps(i,1)=p(l,1)
161 ub(i,1)=ub(1,1)
162 vb(i,1)=vb(1,1)
163 ssou(i,1)=dragu(l)+liftu(l)
164 ssov(i,1)=dragv(l)+liftv(l)
165 blsu(i,1)=phyu(l)-dragu(l)-liftu(l)
166 blsv(i,1)=phyv(l)-dragv(l)-liftv(l)
167
168 enddo
169
170
171 do j = 2,jjm
172
173 C Values at Greenwich (Periodicity)
174
175 zs(iim+1,j)=phis(l+1)/rg
176 ps(iim+1,j)=p(l+1,1)
177 ssou(iim+1,j)=dragu(l+1)+liftu(l+1)
178 ssov(iim+1,j)=dragv(l+1)+liftv(l+1)
179 blsu(iim+1,j)=phyu(l+1)-dragu(l+1)-liftu(l+1)
180 blsv(iim+1,j)=phyv(l+1)-dragv(l+1)-liftv(l+1)
181 zlon(iim+1)=-plon(l+1)*xpi/180.
182 zlat(j)=plat(l+1)*xpi/180.
183
184 ub(iim+1,j)=0.
185 vb(iim+1,j)=0.
186 do k=1,nlev
187 ub(iim+1,j)=ub(iim+1,j)+u(l+1,k)*(p(l+1,k)-p(l+1,k+1))/rg
188 vb(iim+1,j)=vb(iim+1,j)+v(l+1,k)*(p(l+1,k)-p(l+1,k+1))/rg
189 enddo
190
191
192 do i=1,iim
193
194 l=l+1
195 zs(i,j)=phis(l)/rg
196 ps(i,j)=p(l,1)
197 ssou(i,j)=dragu(l)+liftu(l)
198 ssov(i,j)=dragv(l)+liftv(l)
199 blsu(i,j)=phyu(l)-dragu(l)-liftu(l)
200 blsv(i,j)=phyv(l)-dragv(l)-liftv(l)
201 zlon(i)=plon(l)*xpi/180.
202
203 ub(i,j)=0.
204 vb(i,j)=0.
205 do k=1,nlev
206 ub(i,j)=ub(i,j)+u(l,k)*(p(l,k)-p(l,k+1))/rg
207 vb(i,j)=vb(i,j)+v(l,k)*(p(l,k)-p(l,k+1))/rg
208 enddo
209
210 enddo
211
212 enddo
213
214
215 C South Pole
216
217 l=l+1
218 ub(1,jjm+1)=0.
219 vb(1,jjm+1)=0.
220 do k=1,nlev
221 ub(1,jjm+1)=ub(1,jjm+1)+u(l,k)*(p(l,k)-p(l,k+1))/rg
222 vb(1,jjm+1)=vb(1,jjm+1)+v(l,k)*(p(l,k)-p(l,k+1))/rg
223 enddo
224 zlat(jjm+1)=plat(l)*xpi/180.
225
226 do i=1,iim+1
227 zs(i,jjm+1)=phis(l)/rg
228 ps(i,jjm+1)=p(l,1)
229 ssou(i,jjm+1)=dragu(l)+liftu(l)
230 ssov(i,jjm+1)=dragv(l)+liftv(l)
231 blsu(i,jjm+1)=phyu(l)-dragu(l)-liftu(l)
232 blsv(i,jjm+1)=phyv(l)-dragv(l)-liftv(l)
233 ub(i,jjm+1)=ub(1,jjm+1)
234 vb(i,jjm+1)=vb(1,jjm+1)
235 enddo
236
237 C
238 C MOMENT ANGULAIRE
239 C
240 DO j=1,jjm
241 DO i=1,iim
242
243 raam(1)=raam(1)-rea**3*dlon*dlat*0.5*
244 c (cos(zlon(i ))*sin(zlat(j ))*cos(zlat(j ))*ub(i ,j )
245 c +cos(zlon(i ))*sin(zlat(j+1))*cos(zlat(j+1))*ub(i ,j+1))
246 c +rea**3*dlon*dlat*0.5*
247 c (sin(zlon(i ))*cos(zlat(j ))*vb(i ,j )
248 c +sin(zlon(i ))*cos(zlat(j+1))*vb(i ,j+1))
249
250 oaam(1)=oaam(1)-ome*rea**4*dlon*dlat/rg*0.5*
251 c (cos(zlon(i ))*cos(zlat(j ))**2*sin(zlat(j ))*ps(i ,j )
252 c +cos(zlon(i ))*cos(zlat(j+1))**2*sin(zlat(j+1))*ps(i ,j+1))
253
254 raam(2)=raam(2)-rea**3*dlon*dlat*0.5*
255 c (sin(zlon(i ))*sin(zlat(j ))*cos(zlat(j ))*ub(i ,j )
256 c +sin(zlon(i ))*sin(zlat(j+1))*cos(zlat(j+1))*ub(i ,j+1))
257 c -rea**3*dlon*dlat*0.5*
258 c (cos(zlon(i ))*cos(zlat(j ))*vb(i ,j )
259 c +cos(zlon(i ))*cos(zlat(j+1))*vb(i ,j+1))
260
261 oaam(2)=oaam(2)-ome*rea**4*dlon*dlat/rg*0.5*
262 c (sin(zlon(i ))*cos(zlat(j ))**2*sin(zlat(j ))*ps(i ,j )
263 c +sin(zlon(i ))*cos(zlat(j+1))**2*sin(zlat(j+1))*ps(i ,j+1))
264
265 raam(3)=raam(3)+rea**3*dlon*dlat*0.5*
266 c (cos(zlat(j))**2*ub(i,j)+cos(zlat(j+1))**2*ub(i,j+1))
267
268 oaam(3)=oaam(3)+ome*rea**4*dlon*dlat/rg*0.5*
269 c (cos(zlat(j))**3*ps(i,j)+cos(zlat(j+1))**3*ps(i,j+1))
270
271 ENDDO
272 ENDDO
273
274 C
275 C COUPLE DES MONTAGNES:
276 C
277
278 DO j=1,jjm
279 DO i=1,iim
280 tmou(1)=tmou(1)-rea**2*dlon*0.5*sin(zlon(i))
281 c *(zs(i,j)-zs(i,j+1))
282 c *(cos(zlat(j+1))*ps(i,j+1)+cos(zlat(j))*ps(i,j))
283 tmou(2)=tmou(2)+rea**2*dlon*0.5*cos(zlon(i))
284 c *(zs(i,j)-zs(i,j+1))
285 c *(cos(zlat(j+1))*ps(i,j+1)+cos(zlat(j))*ps(i,j))
286 ENDDO
287 ENDDO
288
289 DO j=2,jjm
290 DO i=1,iim
291 tmou(1)=tmou(1)+rea**2*dlat*0.5*sin(zlat(j))
292 c *(zs(i+1,j)-zs(i,j))
293 c *(cos(zlon(i+1))*ps(i+1,j)+cos(zlon(i))*ps(i,j))
294 tmou(2)=tmou(2)+rea**2*dlat*0.5*sin(zlat(j))
295 c *(zs(i+1,j)-zs(i,j))
296 c *(sin(zlon(i+1))*ps(i+1,j)+sin(zlon(i))*ps(i,j))
297 tmou(3)=tmou(3)-rea**2*dlat*0.5*
298 c cos(zlat(j))*(zs(i+1,j)-zs(i,j))*(ps(i+1,j)+ps(i,j))
299 ENDDO
300 ENDDO
301
302 C
303 C COUPLES DES DIFFERENTES FRICTION AU SOL:
304 C
305 l=1
306 DO j=2,jjm
307 DO i=1,iim
308 l=l+1
309 tsso(1)=tsso(1)-rea**3*cos(zlat(j))*dlon*dlat*
310 c ssou(i,j) *sin(zlat(j))*cos(zlon(i))
311 c +rea**3*cos(zlat(j))*dlon*dlat*
312 c ssov(i,j) *sin(zlon(i))
313
314 tsso(2)=tsso(2)-rea**3*cos(zlat(j))*dlon*dlat*
315 c ssou(i,j) *sin(zlat(j))*sin(zlon(i))
316 c -rea**3*cos(zlat(j))*dlon*dlat*
317 c ssov(i,j) *cos(zlon(i))
318
319 tsso(3)=tsso(3)+rea**3*cos(zlat(j))*dlon*dlat*
320 c ssou(i,j) *cos(zlat(j))
321
322 tbls(1)=tbls(1)-rea**3*cos(zlat(j))*dlon*dlat*
323 c blsu(i,j) *sin(zlat(j))*cos(zlon(i))
324 c +rea**3*cos(zlat(j))*dlon*dlat*
325 c blsv(i,j) *sin(zlon(i))
326
327 tbls(2)=tbls(2)-rea**3*cos(zlat(j))*dlon*dlat*
328 c blsu(i,j) *sin(zlat(j))*sin(zlon(i))
329 c -rea**3*cos(zlat(j))*dlon*dlat*
330 c blsv(i,j) *cos(zlon(i))
331
332 tbls(3)=tbls(3)+rea**3*cos(zlat(j))*dlon*dlat*
333 c blsu(i,j) *cos(zlat(j))
334
335 ENDDO
336 ENDDO
337
338
339 100 format(F12.5,15(1x,F12.5))
340
341 aam=raam(3)
342 torsfc= tmou(3)+tsso(3)+tbls(3)
343 c
344 RETURN
345 END

  ViewVC Help
Powered by ViewVC 1.1.21