/[lmdze]/trunk/phylmd/aaam_bud.f90
ViewVC logotype

Diff of /trunk/phylmd/aaam_bud.f90

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/libf/phylmd/aaam_bud.f revision 3 by guez, Wed Feb 27 13:16:39 2008 UTC trunk/libf/phylmd/aaam_bud.f90 revision 61 by guez, Fri Apr 20 14:58:43 2012 UTC
# Line 1  Line 1 
1        subroutine aaam_bud (iam,nlon,nlev,rjour,rsec,  module aaam_bud_m
2       i                   rea,rg,ome,        
3       i                   plat,plon,phis,    implicit none
4       i                   dragu,liftu,phyu,  
5       i                   dragv,liftv,phyv,  contains
6       i                   p, u, v,  
7       o                   aam, torsfc)    subroutine aaam_bud(rea, rg, ome, plat, plon, phis, dragu, liftu, phyu, &
8  c         dragv, liftv, phyv, p, u, v, aam, torsfc)
9        use dimens_m  
10        use dimphy      ! Author: F. Lott (LMD/CNRS). Date: 2003/10/20. Object: Compute
11        implicit none      ! different terms of the axial AAAM budget and mountain torque.
12  c======================================================================      ! Only valid for regular rectangular grids. Should be called after
13  c Auteur(s): F.Lott (LMD/CNRS) date: 20031020      ! "lift_noro".
14  c Object: Compute different terms of the axial AAAM Budget.  
15  C No outputs, every AAM quantities are written on the IAM      USE dimens_m, ONLY : iim, jjm
16  C File.      use nr_util, only: assert_eq, assert, pi
17  c  
18  c Modif : I.Musat (LMD/CNRS) date : 20041020      real, intent(in):: rea ! Earth radius
19  c Outputs : axial components of wind AAM "aam" and total surface torque "torsfc",      real, intent(in):: rg ! gravity constant
20  c but no write in the iam file.      real, intent(in):: ome ! Earth rotation rate
21  c  
22  C WARNING: Only valid for regular rectangular grids.      REAL, intent(in):: plat(:), plon(:)
23  C REMARK: CALL DANS PHYSIQ AFTER lift_noro:      ! (nlon) latitude and longitude in degrees
24  C        CALL aaam_bud (27,klon,klev,rjourvrai,gmtime,  
25  C    C               ra,rg,romega,      real, intent(in):: phis(:) ! (nlon) Geopotential at the ground
26  C    C               rlat,rlon,pphis,      REAL, intent(in):: dragu(:) ! (nlon) orodrag stress (zonal)
27  C    C               zustrdr,zustrli,zustrph,      REAL, intent(in):: liftu(:) ! (nlon) orolift stress (zonal)
28  C    C               zvstrdr,zvstrli,zvstrph,      REAL, intent(in):: phyu(:) ! (nlon) Stress total de la physique (zonal)
29  C    C               paprs,u,v)      REAL, intent(in):: dragv(:) ! (nlon) orodrag stress (Meridional)
30  C      REAL, intent(in):: liftv(:) ! (nlon) orolift stress (Meridional)
31  C======================================================================      REAL, intent(in):: phyv(:) ! (nlon) Stress total de la physique (Meridional)
32  c Explicit Arguments:  
33  c ==================      REAL, intent(in):: p(:, :)
34  c iam-----input-I-File number where AAMs and torques are written      ! (nlon, nlev + 1) pressure (Pa) at model half levels
35  c                 It is a formatted file that has been opened  
36  c                 in physiq.F      real, intent(in):: u(:, :), v(:, :) ! (nlon, nlev) horizontal wind (m/s)
37  c nlon----input-I-Total number of horizontal points that get into physics      REAL, intent(out):: aam ! axial component of wind AAM
38  c nlev----input-I-Number of vertical levels      REAL, intent(out):: torsfc ! axial component of total surface torque
39  c rjour---input-R-Jour compte depuis le debut de la simu (run.def)  
40  c rsec----input-R-Seconde de la journee      ! Local Variables:
41  c rea-----input-R-Earth radius  
42  c rg------input-R-gravity constant      INTEGER nlev ! number of vertical levels
43  c ome-----input-R-Earth rotation rate      INTEGER i, j, k, l
44  c plat ---input-R-Latitude en degres      REAL hadley, hadday
45  c plon ---input-R-Longitude en degres      REAL dlat, dlon ! latitude and longitude increments (radians)
46  c phis ---input-R-Geopotential at the ground  
47  c dragu---input-R-orodrag stress (zonal)      REAL raam(3) ! wind AAM (components 1 & 2: equatorial; component 3: axial)
48  c liftu---input-R-orolift stress (zonal)      REAL oaam(3) ! mass AAM (components 1 & 2: equatorial; component 3: axial)
49  c phyu----input-R-Stress total de la physique (zonal)      REAL tmou(3) ! resolved mountain torque (3 components)
50  c dragv---input-R-orodrag stress (Meridional)      REAL tsso(3) ! parameterised moutain drag torque (3 components)
51  c liftv---input-R-orolift stress (Meridional)      REAL tbls(3) ! parameterised boundary layer torque (3 components)
52  c phyv----input-R-Stress total de la physique (Meridional)      integer iax
53  c p-------input-R-Pressure (Pa) at model half levels  
54  c u-------input-R-Horizontal wind (m/s)      REAL ZS(801, 401) ! topographic height
55  c v-------input-R-Meridional wind (m/s)      REAL PS(801, 401) ! surface pressure
56  c aam-----output-R-Axial Wind AAM (=raam(3))      REAL UB(801, 401), VB(801, 401) ! barotropic wind, zonal and meridional
57  c torsfc--output-R-Total surface torque (=tmou(3)+tsso(3)+tbls(3))      REAL SSOU(801, 401), SSOV(801, 401)
58  c      REAL BLSU(801, 401), BLSV(801, 401)
59  c Implicit Arguments:      REAL ZLON(801), ZLAT(401) ! longitude and latitude in radians
60  c ===================  
61  c      !-------------------------------------------------------------------
62  c iim--common-I: Number of longitude intervals  
63  c jjm--common-I: Number of latitude intervals      call assert(size(plat) == (/size(plon), size(phis), size(dragu), &
64  c klon-common-I: Number of points seen by the physics           size(liftu), size(phyu), size(dragv), size(liftv), size(phyv), &
65  c                iim*(jjm-1)+2 for instance           size(p, 1), size(u, 1), size(v, 1)/), "aaam_bud nlon")
66  c klev-common-I: Number of vertical layers      nlev = assert_eq(size(p, 2) - 1, size(u, 2), size(v, 2), "aaam_bud nlev")
67  c======================================================================  
68  c Local Variables:      if (iim + 1 > 801 .or. jjm + 1 > 401) then
69  c ================         print *, ' Problème de dimension dans aaam_bud'
70  c dlat-----R: Latitude increment (Radians)         stop 1
71  c dlon-----R: Longitude increment (Radians)      endif
72  c raam  ---R: Wind AAM (3 Components, 1 & 2 Equatoriales; 3 Axiale)  
73  c oaam  ---R: Mass AAM (3 Components, 1 & 2 Equatoriales; 3 Axiale)      hadley = 1e18
74  c tmou-----R: Resolved Mountain torque (3 components)      hadday = 1e18 * 24. * 3600.
75  c tsso-----R: Parameterised Moutain drag torque (3 components)      dlat = pi / jjm
76  c tbls-----R: Parameterised Boundary layer torque (3 components)      dlon = 2 * pi / real(iim)
77  c  
78  c LOCAL ARRAY:      oaam = 0.
79  c ===========      raam = 0.
80  c zs    ---R: Topographic height      tmou = 0.
81  c ps    ---R: Surface Pressure        tsso = 0.
82  c ub    ---R: Barotropic wind zonal      tbls = 0.
83  c vb    ---R: Barotropic wind meridional  
84  c zlat  ---R: Latitude in radians      ! Mountain height, pressure and barotropic wind:
85  c zlon  ---R: Longitude in radians  
86  c======================================================================      ! North pole values (j = 1):
87    
88  c      ub(1, 1) = 0.
89  c ARGUMENTS      vb(1, 1) = 0.
90  c      do k = 1, nlev
91        INTEGER iam,nlon,nlev         ub(1, 1) = ub(1, 1) + u(1, k) * (p(1, k) - p(1, k + 1)) / rg
92        REAL rjour         vb(1, 1) = vb(1, 1) + v(1, k) * (p(1, k) - p(1, k + 1)) / rg
93        real, intent(in):: rsec      enddo
94        real rea  
95        real, intent(in):: rg      zlat(1) = plat(1) * pi / 180.
96        real ome  
97        REAL, intent(in):: plat(nlon),plon(nlon)      do i = 1, iim + 1
98        real phis(nlon)         zs(i, 1) = phis(1) / rg
99        REAL dragu(nlon),liftu(nlon),phyu(nlon)                     ps(i, 1) = p(1, 1)
100        REAL dragv(nlon),liftv(nlon),phyv(nlon)                     ub(i, 1) = ub(1, 1)
101        REAL, intent(in):: p(nlon,nlev+1)         vb(i, 1) = vb(1, 1)
102        real u(nlon,nlev), v(nlon,nlev)         ssou(i, 1) = dragu(1) + liftu(1)
103  c         ssov(i, 1) = dragv(1) + liftv(1)
104  c Variables locales:         blsu(i, 1) = phyu(1) - dragu(1) - liftu(1)
105  c         blsv(i, 1) = phyv(1) - dragv(1) - liftv(1)
106        INTEGER i,j,k,l      enddo
107        REAL xpi,hadley,hadday  
108        REAL dlat,dlon      l = 1
109        REAL raam(3),oaam(3),tmou(3),tsso(3),tbls(3)      do j = 2, jjm
110        integer iax         ! Values at Greenwich (Periodicity)
111  cIM ajout aam, torsfc  
112  c aam = composante axiale du Wind AAM raam         zs(iim + 1, j) = phis(l + 1) / rg
113  c torsfc = composante axiale de (tmou+tsso+tbls)         ps(iim + 1, j) = p(l + 1, 1)
114        REAL aam, torsfc         ssou(iim + 1, j) = dragu(l + 1) + liftu(l + 1)
115           ssov(iim + 1, j) = dragv(l + 1) + liftv(l + 1)
116        REAL ZS(801,401),PS(801,401)         blsu(iim + 1, j) = phyu(l + 1) - dragu(l + 1) - liftu(l + 1)
117        REAL UB(801,401),VB(801,401)         blsv(iim + 1, j) = phyv(l + 1) - dragv(l + 1) - liftv(l + 1)
118        REAL SSOU(801,401),SSOV(801,401)         zlon(iim + 1) = - plon(l + 1) * pi / 180.
119        REAL BLSU(801,401),BLSV(801,401)         zlat(j) = plat(l + 1) * pi / 180.
120        REAL ZLON(801),ZLAT(401)  
121  C         ub(iim + 1, j) = 0.
122  C  PUT AAM QUANTITIES AT ZERO:         vb(iim + 1, j) = 0.
123  C         do k = 1, nlev
124        if(iim+1.gt.801.or.jjm+1.gt.401)then            ub(iim + 1, j) = ub(iim + 1, j) &
125        print *,' Pb de dimension dans aaam_bud'                 + u(l + 1, k) * (p(l + 1, k) - p(l + 1, k + 1)) / rg
126        stop            vb(iim + 1, j) = vb(iim + 1, j) &
127        endif                 + v(l + 1, k) * (p(l + 1, k) - p(l + 1, k + 1)) / rg
128           enddo
129        xpi=acos(-1.)  
130        hadley=1.e18         do i = 1, iim
131        hadday=1.e18*24.*3600.            l = l + 1
132        dlat=xpi/float(jjm)            zs(i, j) = phis(l) / rg
133        dlon=2.*xpi/float(iim)            ps(i, j) = p(l, 1)
134                    ssou(i, j) = dragu(l) + liftu(l)
135        do iax=1,3            ssov(i, j) = dragv(l) + liftv(l)
136        oaam(iax)=0.            blsu(i, j) = phyu(l) - dragu(l) - liftu(l)
137        raam(iax)=0.            blsv(i, j) = phyv(l) - dragv(l) - liftv(l)
138        tmou(iax)=0.            zlon(i) = plon(l) * pi / 180.
139        tsso(iax)=0.  
140        tbls(iax)=0.            ub(i, j) = 0.
141        enddo            vb(i, j) = 0.
142              do k = 1, nlev
143  C MOUNTAIN HEIGHT, PRESSURE AND BAROTROPIC WIND:               ub(i, j) = ub(i, j) + u(l, k) * (p(l, k) - p(l, k + 1)) / rg
144                 vb(i, j) = vb(i, j) + v(l, k) * (p(l, k) - p(l, k + 1)) / rg
145  C North pole values (j=1):            enddo
146           enddo
147        l=1      enddo
148    
149          ub(1,1)=0.      ! South Pole
150          vb(1,1)=0.  
151          do k=1,nlev      l = l + 1
152            ub(1,1)=ub(1,1)+u(l,k)*(p(l,k)-p(l,k+1))/rg      ub(1, jjm + 1) = 0.
153            vb(1,1)=vb(1,1)+v(l,k)*(p(l,k)-p(l,k+1))/rg      vb(1, jjm + 1) = 0.
154          enddo      do k = 1, nlev
155           ub(1, jjm + 1) = ub(1, jjm + 1) + u(l, k) * (p(l, k) - p(l, k + 1)) / rg
156            zlat(1)=plat(l)*xpi/180.         vb(1, jjm + 1) = vb(1, jjm + 1) + v(l, k) * (p(l, k) - p(l, k + 1)) / rg
157        enddo
158          do i=1,iim+1      zlat(jjm + 1) = plat(l) * pi / 180.
159    
160            zs(i,1)=phis(l)/rg      do i = 1, iim + 1
161            ps(i,1)=p(l,1)         zs(i, jjm + 1) = phis(l) / rg
162            ub(i,1)=ub(1,1)                                     ps(i, jjm + 1) = p(l, 1)
163            vb(i,1)=vb(1,1)                                     ssou(i, jjm + 1) = dragu(l) + liftu(l)
164            ssou(i,1)=dragu(l)+liftu(l)         ssov(i, jjm + 1) = dragv(l) + liftv(l)
165            ssov(i,1)=dragv(l)+liftv(l)         blsu(i, jjm + 1) = phyu(l) - dragu(l) - liftu(l)
166            blsu(i,1)=phyu(l)-dragu(l)-liftu(l)         blsv(i, jjm + 1) = phyv(l) - dragv(l) - liftv(l)
167            blsv(i,1)=phyv(l)-dragv(l)-liftv(l)         ub(i, jjm + 1) = ub(1, jjm + 1)
168           vb(i, jjm + 1) = vb(1, jjm + 1)
169          enddo      enddo
170    
171        ! Moment angulaire
172        do j = 2,jjm  
173        DO j = 1, jjm
174  C Values at Greenwich (Periodicity)         DO i = 1, iim
175              raam(1) = raam(1) - rea**3 * dlon * dlat * 0.5 * (cos(zlon(i )) &
176        zs(iim+1,j)=phis(l+1)/rg                 * sin(zlat(j )) * cos(zlat(j )) * ub(i , j ) + cos(zlon(i )) &
177        ps(iim+1,j)=p(l+1,1)                 * sin(zlat(j + 1)) * cos(zlat(j + 1)) * ub(i , j + 1)) &
178            ssou(iim+1,j)=dragu(l+1)+liftu(l+1)                 + rea**3 * dlon * dlat * 0.5 * (sin(zlon(i )) * cos(zlat(j )) &
179            ssov(iim+1,j)=dragv(l+1)+liftv(l+1)                 * vb(i , j ) + sin(zlon(i )) * cos(zlat(j + 1)) * vb(i , j + 1))
180            blsu(iim+1,j)=phyu(l+1)-dragu(l+1)-liftu(l+1)  
181            blsv(iim+1,j)=phyv(l+1)-dragv(l+1)-liftv(l+1)            oaam(1) = oaam(1) - ome * rea**4 * dlon * dlat / rg * 0.5 &
182        zlon(iim+1)=-plon(l+1)*xpi/180.                 * (cos(zlon(i )) * cos(zlat(j ))**2 * sin(zlat(j )) &
183        zlat(j)=plat(l+1)*xpi/180.                 * ps(i , j ) + cos(zlon(i )) * cos(zlat(j + 1))**2 &
184                   * sin(zlat(j + 1)) * ps(i , j + 1))
185        ub(iim+1,j)=0.  
186        vb(iim+1,j)=0.            raam(2) = raam(2) - rea**3 * dlon * dlat * 0.5 * (sin(zlon(i )) &
187           do k=1,nlev                 * sin(zlat(j )) * cos(zlat(j )) * ub(i , j ) + sin(zlon(i )) &
188           ub(iim+1,j)=ub(iim+1,j)+u(l+1,k)*(p(l+1,k)-p(l+1,k+1))/rg                 * sin(zlat(j + 1)) * cos(zlat(j + 1)) * ub(i , j + 1)) &
189           vb(iim+1,j)=vb(iim+1,j)+v(l+1,k)*(p(l+1,k)-p(l+1,k+1))/rg                 - rea**3 * dlon * dlat * 0.5 * (cos(zlon(i )) * cos(zlat(j )) &
190           enddo                 * vb(i , j ) + cos(zlon(i )) * cos(zlat(j + 1)) * vb(i , j + 1))
191          
192              oaam(2) = oaam(2) - ome * rea**4 * dlon * dlat / rg * 0.5 &
193        do i=1,iim                 * (sin(zlon(i )) * cos(zlat(j ))**2 * sin(zlat(j )) &
194                   * ps(i , j ) + sin(zlon(i )) * cos(zlat(j + 1))**2 &
195        l=l+1                 * sin(zlat(j + 1)) * ps(i , j + 1))
196        zs(i,j)=phis(l)/rg  
197        ps(i,j)=p(l,1)            raam(3) = raam(3) + rea**3 * dlon * dlat * 0.5 * (cos(zlat(j))**2 &
198            ssou(i,j)=dragu(l)+liftu(l)                 * ub(i, j) + cos(zlat(j + 1))**2 * ub(i, j + 1))
199            ssov(i,j)=dragv(l)+liftv(l)  
200            blsu(i,j)=phyu(l)-dragu(l)-liftu(l)            oaam(3) = oaam(3) + ome * rea**4 * dlon * dlat / rg * 0.5 &
201            blsv(i,j)=phyv(l)-dragv(l)-liftv(l)                 * (cos(zlat(j))**3 * ps(i, j) + cos(zlat(j + 1))**3 &
202        zlon(i)=plon(l)*xpi/180.                 * ps(i, j + 1))
203           ENDDO
204        ub(i,j)=0.      ENDDO
205        vb(i,j)=0.  
206           do k=1,nlev      ! Couple des montagnes :
207           ub(i,j)=ub(i,j)+u(l,k)*(p(l,k)-p(l,k+1))/rg  
208           vb(i,j)=vb(i,j)+v(l,k)*(p(l,k)-p(l,k+1))/rg      DO j = 1, jjm
209           enddo         DO i = 1, iim
210              tmou(1) = tmou(1) - rea**2 * dlon * 0.5 * sin(zlon(i)) &
211        enddo                 * (zs(i, j) - zs(i, j + 1)) &
212                   * (cos(zlat(j + 1)) * ps(i, j + 1) + cos(zlat(j)) * ps(i, j))
213        enddo            tmou(2) = tmou(2) + rea**2 * dlon * 0.5 * cos(zlon(i)) &
214                   * (zs(i, j) - zs(i, j + 1)) &
215                   * (cos(zlat(j + 1)) * ps(i, j + 1) + cos(zlat(j)) * ps(i, j))
216  C South Pole         ENDDO
217        ENDDO
218        l=l+1  
219        ub(1,jjm+1)=0.      DO j = 2, jjm
220        vb(1,jjm+1)=0.         DO i = 1, iim
221        do k=1,nlev            tmou(1) = tmou(1) + rea**2 * dlat * 0.5 * sin(zlat(j)) &
222           ub(1,jjm+1)=ub(1,jjm+1)+u(l,k)*(p(l,k)-p(l,k+1))/rg                 * (zs(i + 1, j) - zs(i, j)) &
223           vb(1,jjm+1)=vb(1,jjm+1)+v(l,k)*(p(l,k)-p(l,k+1))/rg                 * (cos(zlon(i + 1)) * ps(i + 1, j) + cos(zlon(i)) * ps(i, j))
224        enddo            tmou(2) = tmou(2) + rea**2 * dlat * 0.5 * sin(zlat(j)) &
225        zlat(jjm+1)=plat(l)*xpi/180.                 * (zs(i + 1, j) - zs(i, j)) &
226                   * (sin(zlon(i + 1)) * ps(i + 1, j) + sin(zlon(i)) * ps(i, j))
227        do i=1,iim+1            tmou(3) = tmou(3) - rea**2 * dlat * 0.5* cos(zlat(j)) &
228        zs(i,jjm+1)=phis(l)/rg                 * (zs(i + 1, j) - zs(i, j)) * (ps(i + 1, j) + ps(i, j))
229        ps(i,jjm+1)=p(l,1)         ENDDO
230            ssou(i,jjm+1)=dragu(l)+liftu(l)      ENDDO
231            ssov(i,jjm+1)=dragv(l)+liftv(l)  
232            blsu(i,jjm+1)=phyu(l)-dragu(l)-liftu(l)      ! Couples des differentes friction au sol :
233            blsv(i,jjm+1)=phyv(l)-dragv(l)-liftv(l)  
234        ub(i,jjm+1)=ub(1,jjm+1)                                    DO j = 2, jjm
235        vb(i,jjm+1)=vb(1,jjm+1)                                         DO i = 1, iim
236        enddo            tsso(1) = tsso(1) - rea**3 * cos(zlat(j)) * dlon * dlat* &
237                   ssou(i, j) * sin(zlat(j)) * cos(zlon(i)) &
238  C                 + rea**3 * cos(zlat(j)) * dlon * dlat* &
239  C  MOMENT ANGULAIRE                 ssov(i, j) * sin(zlon(i))
240  C  
241          DO j=1,jjm                tsso(2) = tsso(2) - rea**3 * cos(zlat(j)) * dlon * dlat* &
242          DO i=1,iim                 ssou(i, j) * sin(zlat(j)) * sin(zlon(i)) &
243                   - rea**3 * cos(zlat(j)) * dlon * dlat* &
244             raam(1)=raam(1)-rea**3*dlon*dlat*0.5*                 ssov(i, j) * cos(zlon(i))
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))            tsso(3) = tsso(3) + rea**3 * cos(zlat(j)) * dlon * dlat* &
247       c                    +rea**3*dlon*dlat*0.5*                 ssou(i, j) * cos(zlat(j))
248       c    (sin(zlon(i  ))*cos(zlat(j  ))*vb(i  ,j  )  
249       c    +sin(zlon(i  ))*cos(zlat(j+1))*vb(i  ,j+1))            tbls(1) = tbls(1) - rea**3 * cos(zlat(j)) * dlon * dlat* &
250                   blsu(i, j) * sin(zlat(j)) * cos(zlon(i)) &
251             oaam(1)=oaam(1)-ome*rea**4*dlon*dlat/rg*0.5*                 + rea**3 * cos(zlat(j)) * dlon * dlat* &
252       c   (cos(zlon(i  ))*cos(zlat(j  ))**2*sin(zlat(j  ))*ps(i  ,j  )                 blsv(i, j) * sin(zlon(i))
253       c   +cos(zlon(i  ))*cos(zlat(j+1))**2*sin(zlat(j+1))*ps(i  ,j+1))  
254              tbls(2) = tbls(2) - rea**3 * cos(zlat(j)) * dlon * dlat* &
255             raam(2)=raam(2)-rea**3*dlon*dlat*0.5*                 blsu(i, j) * sin(zlat(j)) * sin(zlon(i)) &
256       c    (sin(zlon(i  ))*sin(zlat(j  ))*cos(zlat(j  ))*ub(i  ,j  )                 - rea**3 * cos(zlat(j)) * dlon * dlat* &
257       c    +sin(zlon(i  ))*sin(zlat(j+1))*cos(zlat(j+1))*ub(i  ,j+1))                 blsv(i, j) * cos(zlon(i))
258       c                    -rea**3*dlon*dlat*0.5*  
259       c    (cos(zlon(i  ))*cos(zlat(j  ))*vb(i  ,j  )            tbls(3) = tbls(3) + rea**3 * cos(zlat(j)) * dlon * dlat* &
260       c    +cos(zlon(i  ))*cos(zlat(j+1))*vb(i  ,j+1))                 blsu(i, j) * cos(zlat(j))
261           ENDDO
262             oaam(2)=oaam(2)-ome*rea**4*dlon*dlat/rg*0.5*      ENDDO
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))      aam = raam(3)
265        torsfc = tmou(3) + tsso(3) + tbls(3)
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))    END subroutine aaam_bud
268    
269             oaam(3)=oaam(3)+ome*rea**4*dlon*dlat/rg*0.5*  end module aaam_bud_m
      c        (cos(zlat(j))**3*ps(i,j)+cos(zlat(j+1))**3*ps(i,j+1))  
   
         ENDDO  
         ENDDO  
   
 C  
 C COUPLE DES MONTAGNES:  
 C  
   
         DO j=1,jjm  
         DO i=1,iim  
            tmou(1)=tmou(1)-rea**2*dlon*0.5*sin(zlon(i))  
      c  *(zs(i,j)-zs(i,j+1))  
      c  *(cos(zlat(j+1))*ps(i,j+1)+cos(zlat(j))*ps(i,j))  
            tmou(2)=tmou(2)+rea**2*dlon*0.5*cos(zlon(i))  
      c  *(zs(i,j)-zs(i,j+1))  
      c  *(cos(zlat(j+1))*ps(i,j+1)+cos(zlat(j))*ps(i,j))  
         ENDDO  
         ENDDO  
             
         DO j=2,jjm  
         DO i=1,iim  
            tmou(1)=tmou(1)+rea**2*dlat*0.5*sin(zlat(j))  
      c  *(zs(i+1,j)-zs(i,j))  
      c  *(cos(zlon(i+1))*ps(i+1,j)+cos(zlon(i))*ps(i,j))  
            tmou(2)=tmou(2)+rea**2*dlat*0.5*sin(zlat(j))  
      c  *(zs(i+1,j)-zs(i,j))  
      c  *(sin(zlon(i+1))*ps(i+1,j)+sin(zlon(i))*ps(i,j))  
            tmou(3)=tmou(3)-rea**2*dlat*0.5*  
      c  cos(zlat(j))*(zs(i+1,j)-zs(i,j))*(ps(i+1,j)+ps(i,j))  
         ENDDO  
         ENDDO  
   
 C  
 C COUPLES DES DIFFERENTES FRICTION AU SOL:  
 C  
         l=1  
         DO j=2,jjm  
         DO i=1,iim  
         l=l+1  
            tsso(1)=tsso(1)-rea**3*cos(zlat(j))*dlon*dlat*  
      c     ssou(i,j)          *sin(zlat(j))*cos(zlon(i))  
      c                    +rea**3*cos(zlat(j))*dlon*dlat*  
      c     ssov(i,j)          *sin(zlon(i))  
   
            tsso(2)=tsso(2)-rea**3*cos(zlat(j))*dlon*dlat*  
      c     ssou(i,j)          *sin(zlat(j))*sin(zlon(i))  
      c                    -rea**3*cos(zlat(j))*dlon*dlat*  
      c     ssov(i,j)          *cos(zlon(i))  
   
            tsso(3)=tsso(3)+rea**3*cos(zlat(j))*dlon*dlat*  
      c     ssou(i,j)          *cos(zlat(j))  
   
            tbls(1)=tbls(1)-rea**3*cos(zlat(j))*dlon*dlat*  
      c     blsu(i,j)          *sin(zlat(j))*cos(zlon(i))  
      c                    +rea**3*cos(zlat(j))*dlon*dlat*  
      c     blsv(i,j)          *sin(zlon(i))  
   
            tbls(2)=tbls(2)-rea**3*cos(zlat(j))*dlon*dlat*  
      c     blsu(i,j)          *sin(zlat(j))*sin(zlon(i))  
      c                    -rea**3*cos(zlat(j))*dlon*dlat*  
      c     blsv(i,j)          *cos(zlon(i))  
   
            tbls(3)=tbls(3)+rea**3*cos(zlat(j))*dlon*dlat*  
      c     blsu(i,j)          *cos(zlat(j))  
   
         ENDDO  
         ENDDO  
               
   
 100   format(F12.5,15(1x,F12.5))  
   
       aam=raam(3)  
       torsfc= tmou(3)+tsso(3)+tbls(3)  
 c  
       RETURN  
       END  

Legend:
Removed from v.3  
changed lines
  Added in v.61

  ViewVC Help
Powered by ViewVC 1.1.21