source: trunk/LMDZ.TITAN/libf/phytitan/aaam_bud.F @ 1243

Last change on this file since 1243 was 102, checked in by slebonnois, 14 years ago

SL : corrections et modifications dans phytitan correspondant a celles
faites apres compilation Venus. Titan pas encore compile.

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      use dimphy
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 WARNING: Only valid for regular rectangular grids.
16C REMARK: CALL DANS PHYSIQ AFTER lift_noro:
17C        CALL aaam_bud (27,klon,klev,rjourvrai,gmtime,
18C    C               ra,rg,romega,
19C    C               rlat,rlon,pphis,
20C    C               zustrdr,zustrli,zustrcl,
21C    C               zvstrdr,zvstrli,zvstrcl,
22C    C               paprs,u,v)
23C
24C======================================================================
25c Explicit Arguments:
26c ==================
27c iam-----input-I-File number where AAMs and torques are written
28c                 It is a formatted file that has been opened
29c                 in physiq.F
30c nlon----input-I-Total number of horizontal points that get into physics
31c nlev----input-I-Number of vertical levels
32c rjour---input-R-Jour compte depuis le debut de la simu (run.def)
33c rsec----input-R-Seconde de la journee
34c rea-----input-R-Earth radius
35c rg------input-R-gravity constant
36c ome-----input-R-Earth rotation rate
37c plat ---input-R-Latitude en degres
38c plon ---input-R-Longitude en degres
39c phis ---input-R-Geopotential at the ground
40c dragu---input-R-orodrag stress (zonal)
41c liftu---input-R-orolift stress (zonal)
42c clu-----input-R-Boundary layer stress (zonal)
43c dragv---input-R-orodrag stress (Meridional)
44c liftv---input-R-orolift stress (Meridional)
45c clv-----input-R-Boundary layer stress (Meridional)
46c p-------input-R-Pressure (Pa) at model half levels
47c u-------input-R-Horizontal wind (m/s)
48c v-------input-R-Meridional wind (m/s)
49c
50c
51c Implicit Arguments:
52c ===================
53c
54c iim--common-I: Number of longitude intervals
55c jjm--common-I: Number of latitude intervals
56c klon-common-I: Number of points seen by the physics
57c                iim*(jjm-1)+2 for instance
58c klev-common-I: Number of vertical layers
59c======================================================================
60c Local Variables:
61c ================
62c dlat-----R: Latitude increment (Radians)
63c dlon-----R: Longitude increment (Radians)
64c raam  ---R: Wind AAM (3 Components, 1 & 2 Equatoriales; 3 Axiale)
65c oaam  ---R: Mass AAM (3 Components, 1 & 2 Equatoriales; 3 Axiale)
66c tmou-----R: Resolved Mountain torque (3 components)
67c tsso-----R: Parameterised Moutain drag torque (3 components)
68c tbls-----R: Parameterised Boundary layer torque (3 components)
69c
70c LOCAL ARRAY:
71c ===========
72c zs    ---R: Topographic height
73c ps    ---R: Surface Pressure 
74c ub    ---R: Barotropic wind zonal
75c vb    ---R: Barotropic wind meridional
76c zlat  ---R: Latitude in radians
77c zlon  ---R: Longitude in radians
78c======================================================================
79
80#include "dimensions.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,
328     c      raam(1)/hadday,oaam(1)/hadday,
329     c      tmou(1)/hadley,tsso(1)/hadley,tbls(1)/hadley,
330     c      raam(2)/hadday,oaam(2)/hadday,
331     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
334100   format(F12.5,15(1x,F12.5))
335c00   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.