source: trunk/libf/phyvenus/aaam_bud.F @ 24

Last change on this file since 24 was 3, checked in by slebonnois, 14 years ago

Creation de repertoires:

  • chantiers : pour communiquer sur nos projets de modifs
  • documentation : pour stocker les docs

Ajout de:

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