source: LMDZ4/trunk/libf/phy_IPCC_AR4/aaam_bud.F @ 941

Last change on this file since 941 was 868, checked in by Laurent Fairhead, 17 years ago

Preparation du remplacement de la physique utilisee pour l'exercice IPCC_AR4
par la version de la physique avec thermique. On garde le repertoire phylmd
pour un petit moment pour que les utilisateurs ne soient pas trop perdus ...
phy_IPCC_AR4 = phylmd
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.