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

Last change on this file since 1357 was 1279, checked in by Laurent Fairhead, 15 years ago

Merged LMDZ4-dev branch changes r1241:1278 into the trunk
Running trunk and LMDZ4-dev in LMDZOR configuration on local
machine (sequential) and SX8 (4-proc) yields identical results
(restart and restartphy are identical binarily)
Log history from r1241 to r1278 is available by switching to
source:LMDZ4/branches/LMDZ4-dev-20091210

  • 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        -R-Jour compte depuis le debut de la simu (run.def)
42c rsec         -R-Seconde de la journee
43c rea          -R-Earth radius
44c rg           -R-gravity constant
45c ome          -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, 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)
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.