source: LMDZ4/trunk/libf/phylmd/aaam_bud.F @ 862

Last change on this file since 862 was 644, checked in by Laurent Fairhead, 19 years ago

Synchronisation avec tous les diagnostiques de Ionela IM
Inclusion du slab ocean IM
LF

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