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

Last change on this file since 3303 was 1530, checked in by emillour, 9 years ago

Venus and Titan GCMs:
Updates in the physics to keep up with updates in LMDZ5 (up to
LMDZ5 trunk, rev 2350) concerning dynamics/physics separation:

  • Adapted makelmdz and makelmdz_fcm script to stop if trying to compile 1d model or newstart or start2archive in parallel.
  • got rid of references to "dimensions.h" in physics. Within physics packages, use nbp_lon (=iim), nbp_lat (=jjmp1) and nbp_lev (=llm) from module mod_grid_phy_lmdz (in phy_common) instead. Only partially done for Titan, because of many hard-coded commons; a necessary first step will be to clean these up (using modules).

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