source: LMDZ4/branches/LMDZ4-dev/libf/phylmd/aaam_bud.F @ 1105

Last change on this file since 1105 was 940, checked in by Laurent Fairhead, 16 years ago

On remplace le fichier include dimphy.h par le module dimphy.F90i pour etre
coherent avec le partout
LF

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