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