Sources/phylmd/aaam_bud.f

module aaam_bud_m

  implicit none

contains

  subroutine aaam_bud(rg, ome, plat, plon, phis, dragu, liftu, phyu, dragv, &
       liftv, phyv, p, u, v, aam, torsfc)

    ! Author: F. Lott (LMD/CNRS). Date: 2003/10/20. Object: Compute
    ! different terms of the axial AAAM budget and mountain torque.
    ! Only valid for regular rectangular grids. Should be called after
    ! "lift_noro".

    USE dimens_m, ONLY : iim, jjm
    use nr_util, only: assert_eq, assert, pi
    USE suphec_m, ONLY: ra

    real, intent(in):: rg ! gravity constant
    real, intent(in):: ome ! Earth rotation rate

    REAL, intent(in):: plat(:), plon(:)
    ! (nlon) latitude and longitude in degrees

    real, intent(in):: phis(:) ! (nlon) Geopotential at the ground
    REAL, intent(in):: dragu(:) ! (nlon) orodrag stress (zonal)
    REAL, intent(in):: liftu(:) ! (nlon) orolift stress (zonal)
    REAL, intent(in):: phyu(:) ! (nlon) Stress total de la physique (zonal)
    REAL, intent(in):: dragv(:) ! (nlon) orodrag stress (Meridional)
    REAL, intent(in):: liftv(:) ! (nlon) orolift stress (Meridional)
    REAL, intent(in):: phyv(:) ! (nlon) Stress total de la physique (Meridional)

    REAL, intent(in):: p(:, :) 
    ! (nlon, nlev + 1) pressure (Pa) at model half levels

    real, intent(in):: u(:, :), v(:, :) ! (nlon, nlev) horizontal wind (m/s)
    REAL, intent(out):: aam ! axial component of wind AAM
    REAL, intent(out):: torsfc ! axial component of total surface torque

    ! Local Variables:

    INTEGER nlev ! number of vertical levels
    INTEGER i, j, k, l
    REAL dlat, dlon ! latitude and longitude increments (radians)

    REAL raam(3) ! wind AAM (components 1 & 2: equatorial; component 3: axial)
    REAL oaam(3) ! mass AAM (components 1 & 2: equatorial; component 3: axial)
    REAL tmou(3) ! resolved mountain torque (3 components)
    REAL tsso(3) ! parameterised moutain drag torque (3 components)
    REAL tbls(3) ! parameterised boundary layer torque (3 components)

    REAL ZS(801, 401) ! topographic height
    REAL PS(801, 401) ! surface pressure 
    REAL UB(801, 401), VB(801, 401) ! barotropic wind, zonal and meridional
    REAL SSOU(801, 401), SSOV(801, 401)
    REAL BLSU(801, 401), BLSV(801, 401)
    REAL ZLON(801), ZLAT(401) ! longitude and latitude in radians

    !-------------------------------------------------------------------

    call assert(size(plat) == (/size(plon), size(phis), size(dragu), &
         size(liftu), size(phyu), size(dragv), size(liftv), size(phyv), &
         size(p, 1), size(u, 1), size(v, 1)/), "aaam_bud nlon")
    nlev = assert_eq(size(p, 2) - 1, size(u, 2), size(v, 2), "aaam_bud nlev")

    if (iim + 1 > 801 .or. jjm + 1 > 401) then
       print *, ' Problème de dimension dans aaam_bud'
       stop 1
    endif

    dlat = pi / jjm
    dlon = 2 * pi / real(iim) 

    oaam = 0.
    raam = 0.
    tmou = 0.
    tsso = 0.
    tbls = 0.

    ! Mountain height, pressure and barotropic wind:

    ! North pole values (j = 1):

    ub(1, 1) = 0.
    vb(1, 1) = 0.
    do k = 1, nlev
       ub(1, 1) = ub(1, 1) + u(1, k) * (p(1, k) - p(1, k + 1)) / rg
       vb(1, 1) = vb(1, 1) + v(1, k) * (p(1, k) - p(1, k + 1)) / rg
    enddo

    zlat(1) = plat(1) * pi / 180.

    do i = 1, iim + 1
       zs(i, 1) = phis(1) / rg
       ps(i, 1) = p(1, 1)
       ub(i, 1) = ub(1, 1) 
       vb(i, 1) = vb(1, 1) 
       ssou(i, 1) = dragu(1) + liftu(1)
       ssov(i, 1) = dragv(1) + liftv(1)
       blsu(i, 1) = phyu(1) - dragu(1) - liftu(1)
       blsv(i, 1) = phyv(1) - dragv(1) - liftv(1)
    enddo

    l = 1
    do j = 2, jjm
       ! Values at Greenwich (Periodicity)

       zs(iim + 1, j) = phis(l + 1) / rg
       ps(iim + 1, j) = p(l + 1, 1)
       ssou(iim + 1, j) = dragu(l + 1) + liftu(l + 1)
       ssov(iim + 1, j) = dragv(l + 1) + liftv(l + 1)
       blsu(iim + 1, j) = phyu(l + 1) - dragu(l + 1) - liftu(l + 1)
       blsv(iim + 1, j) = phyv(l + 1) - dragv(l + 1) - liftv(l + 1)
       zlon(iim + 1) = - plon(l + 1) * pi / 180.
       zlat(j) = plat(l + 1) * pi / 180.

       ub(iim + 1, j) = 0.
       vb(iim + 1, j) = 0.
       do k = 1, nlev
          ub(iim + 1, j) = ub(iim + 1, j) &
               + u(l + 1, k) * (p(l + 1, k) - p(l + 1, k + 1)) / rg
          vb(iim + 1, j) = vb(iim + 1, j) &
               + v(l + 1, k) * (p(l + 1, k) - p(l + 1, k + 1)) / rg
       enddo

       do i = 1, iim
          l = l + 1
          zs(i, j) = phis(l) / rg
          ps(i, j) = p(l, 1)
          ssou(i, j) = dragu(l) + liftu(l)
          ssov(i, j) = dragv(l) + liftv(l)
          blsu(i, j) = phyu(l) - dragu(l) - liftu(l)
          blsv(i, j) = phyv(l) - dragv(l) - liftv(l)
          zlon(i) = plon(l) * pi / 180.

          ub(i, j) = 0.
          vb(i, j) = 0.
          do k = 1, nlev
             ub(i, j) = ub(i, j) + u(l, k) * (p(l, k) - p(l, k + 1)) / rg
             vb(i, j) = vb(i, j) + v(l, k) * (p(l, k) - p(l, k + 1)) / rg
          enddo
       enddo
    enddo

    ! South Pole

    l = l + 1
    ub(1, jjm + 1) = 0.
    vb(1, jjm + 1) = 0.
    do k = 1, nlev
       ub(1, jjm + 1) = ub(1, jjm + 1) + u(l, k) * (p(l, k) - p(l, k + 1)) / rg
       vb(1, jjm + 1) = vb(1, jjm + 1) + v(l, k) * (p(l, k) - p(l, k + 1)) / rg
    enddo
    zlat(jjm + 1) = plat(l) * pi / 180.

    do i = 1, iim + 1
       zs(i, jjm + 1) = phis(l) / rg
       ps(i, jjm + 1) = p(l, 1)
       ssou(i, jjm + 1) = dragu(l) + liftu(l)
       ssov(i, jjm + 1) = dragv(l) + liftv(l)
       blsu(i, jjm + 1) = phyu(l) - dragu(l) - liftu(l)
       blsv(i, jjm + 1) = phyv(l) - dragv(l) - liftv(l)
       ub(i, jjm + 1) = ub(1, jjm + 1) 
       vb(i, jjm + 1) = vb(1, jjm + 1) 
    enddo

    ! Moment angulaire 

    DO j = 1, jjm 
       DO i = 1, iim
          raam(1) = raam(1) - ra**3 * dlon * dlat * 0.5 * (cos(zlon(i )) &
               * sin(zlat(j )) * cos(zlat(j )) * ub(i , j ) + cos(zlon(i )) &
               * sin(zlat(j + 1)) * cos(zlat(j + 1)) * ub(i , j + 1)) &
               + ra**3 * dlon * dlat * 0.5 * (sin(zlon(i )) * cos(zlat(j )) &
               * vb(i , j ) + sin(zlon(i )) * cos(zlat(j + 1)) * vb(i , j + 1))

          oaam(1) = oaam(1) - ome * ra**4 * dlon * dlat / rg * 0.5 &
               * (cos(zlon(i )) * cos(zlat(j ))**2 * sin(zlat(j )) &
               * ps(i , j ) + cos(zlon(i )) * cos(zlat(j + 1))**2 &
               * sin(zlat(j + 1)) * ps(i , j + 1))

          raam(2) = raam(2) - ra**3 * dlon * dlat * 0.5 * (sin(zlon(i )) &
               * sin(zlat(j )) * cos(zlat(j )) * ub(i , j ) + sin(zlon(i )) &
               * sin(zlat(j + 1)) * cos(zlat(j + 1)) * ub(i , j + 1)) &
               - ra**3 * dlon * dlat * 0.5 * (cos(zlon(i )) * cos(zlat(j )) &
               * vb(i , j ) + cos(zlon(i )) * cos(zlat(j + 1)) * vb(i , j + 1))

          oaam(2) = oaam(2) - ome * ra**4 * dlon * dlat / rg * 0.5 &
               * (sin(zlon(i )) * cos(zlat(j ))**2 * sin(zlat(j )) &
               * ps(i , j ) + sin(zlon(i )) * cos(zlat(j + 1))**2 &
               * sin(zlat(j + 1)) * ps(i , j + 1))

          raam(3) = raam(3) + ra**3 * dlon * dlat * 0.5 * (cos(zlat(j))**2 &
               * ub(i, j) + cos(zlat(j + 1))**2 * ub(i, j + 1))

          oaam(3) = oaam(3) + ome * ra**4 * dlon * dlat / rg * 0.5 &
               * (cos(zlat(j))**3 * ps(i, j) + cos(zlat(j + 1))**3 &
               * ps(i, j + 1))
       ENDDO
    ENDDO

    ! Couple des montagnes :

    DO j = 1, jjm
       DO i = 1, iim
          tmou(1) = tmou(1) - ra**2 * dlon * 0.5 * sin(zlon(i)) &
               * (zs(i, j) - zs(i, j + 1)) &
               * (cos(zlat(j + 1)) * ps(i, j + 1) + cos(zlat(j)) * ps(i, j)) 
          tmou(2) = tmou(2) + ra**2 * dlon * 0.5 * cos(zlon(i)) &
               * (zs(i, j) - zs(i, j + 1)) &
               * (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) + ra**2 * dlat * 0.5 * sin(zlat(j)) &
               * (zs(i + 1, j) - zs(i, j)) &
               * (cos(zlon(i + 1)) * ps(i + 1, j) + cos(zlon(i)) * ps(i, j))
          tmou(2) = tmou(2) + ra**2 * dlat * 0.5 * sin(zlat(j)) &
               * (zs(i + 1, j) - zs(i, j)) &
               * (sin(zlon(i + 1)) * ps(i + 1, j) + sin(zlon(i)) * ps(i, j))
          tmou(3) = tmou(3) - ra**2 * dlat * 0.5* cos(zlat(j)) &
               * (zs(i + 1, j) - zs(i, j)) * (ps(i + 1, j) + ps(i, j))
       ENDDO
    ENDDO

    ! Couples des differentes friction au sol :

    DO j = 2, jjm
       DO i = 1, iim
          tsso(1) = tsso(1) - ra**3 * cos(zlat(j)) * dlon * dlat* &
               ssou(i, j) * sin(zlat(j)) * cos(zlon(i)) &
               + ra**3 * cos(zlat(j)) * dlon * dlat* &
               ssov(i, j) * sin(zlon(i))

          tsso(2) = tsso(2) - ra**3 * cos(zlat(j)) * dlon * dlat* &
               ssou(i, j) * sin(zlat(j)) * sin(zlon(i)) &
               - ra**3 * cos(zlat(j)) * dlon * dlat* &
               ssov(i, j) * cos(zlon(i))

          tsso(3) = tsso(3) + ra**3 * cos(zlat(j)) * dlon * dlat* &
               ssou(i, j) * cos(zlat(j))

          tbls(1) = tbls(1) - ra**3 * cos(zlat(j)) * dlon * dlat* &
               blsu(i, j) * sin(zlat(j)) * cos(zlon(i)) &
               + ra**3 * cos(zlat(j)) * dlon * dlat* &
               blsv(i, j) * sin(zlon(i))

          tbls(2) = tbls(2) - ra**3 * cos(zlat(j)) * dlon * dlat* &
               blsu(i, j) * sin(zlat(j)) * sin(zlon(i)) &
               - ra**3 * cos(zlat(j)) * dlon * dlat* &
               blsv(i, j) * cos(zlon(i))

          tbls(3) = tbls(3) + ra**3 * cos(zlat(j)) * dlon * dlat* &
               blsu(i, j) * cos(zlat(j))
       ENDDO
    ENDDO

    aam = raam(3)
    torsfc = tmou(3) + tsso(3) + tbls(3)

  END subroutine aaam_bud

end module aaam_bud_m
1	module aaam_bud_m
2
3	implicit none
4
5	contains
6
7	subroutine aaam_bud(rg, ome, plat, plon, phis, dragu, liftu, phyu, dragv, &
8	liftv, phyv, p, u, v, aam, torsfc)
9
10	! Author: F. Lott (LMD/CNRS). Date: 2003/10/20. Object: Compute
11	! different terms of the axial AAAM budget and mountain torque.
12	! Only valid for regular rectangular grids. Should be called after
13	! "lift_noro".
14
15	USE dimens_m, ONLY : iim, jjm
16	use nr_util, only: assert_eq, assert, pi
17	USE suphec_m, ONLY: ra
18
19	real, intent(in):: rg ! gravity constant
20	real, intent(in):: ome ! Earth rotation rate
21
22	REAL, intent(in):: plat(:), plon(:)
23	! (nlon) latitude and longitude in degrees
24
25	real, intent(in):: phis(:) ! (nlon) Geopotential at the ground
26	REAL, intent(in):: dragu(:) ! (nlon) orodrag stress (zonal)
27	REAL, intent(in):: liftu(:) ! (nlon) orolift stress (zonal)
28	REAL, intent(in):: phyu(:) ! (nlon) Stress total de la physique (zonal)
29	REAL, intent(in):: dragv(:) ! (nlon) orodrag stress (Meridional)
30	REAL, intent(in):: liftv(:) ! (nlon) orolift stress (Meridional)
31	REAL, intent(in):: phyv(:) ! (nlon) Stress total de la physique (Meridional)
32
33	REAL, intent(in):: p(:, :)
34	! (nlon, nlev + 1) pressure (Pa) at model half levels
35
36	real, intent(in):: u(:, :), v(:, :) ! (nlon, nlev) horizontal wind (m/s)
37	REAL, intent(out):: aam ! axial component of wind AAM
38	REAL, intent(out):: torsfc ! axial component of total surface torque
39
40	! Local Variables:
41
42	INTEGER nlev ! number of vertical levels
43	INTEGER i, j, k, l
44	REAL dlat, dlon ! latitude and longitude increments (radians)
45
46	REAL raam(3) ! wind AAM (components 1 & 2: equatorial; component 3: axial)
47	REAL oaam(3) ! mass AAM (components 1 & 2: equatorial; component 3: axial)
48	REAL tmou(3) ! resolved mountain torque (3 components)
49	REAL tsso(3) ! parameterised moutain drag torque (3 components)
50	REAL tbls(3) ! parameterised boundary layer torque (3 components)
51
52	REAL ZS(801, 401) ! topographic height
53	REAL PS(801, 401) ! surface pressure
54	REAL UB(801, 401), VB(801, 401) ! barotropic wind, zonal and meridional
55	REAL SSOU(801, 401), SSOV(801, 401)
56	REAL BLSU(801, 401), BLSV(801, 401)
57	REAL ZLON(801), ZLAT(401) ! longitude and latitude in radians
58
59	!-------------------------------------------------------------------
60
61	call assert(size(plat) == (/size(plon), size(phis), size(dragu), &
62	size(liftu), size(phyu), size(dragv), size(liftv), size(phyv), &
63	size(p, 1), size(u, 1), size(v, 1)/), "aaam_bud nlon")
64	nlev = assert_eq(size(p, 2) - 1, size(u, 2), size(v, 2), "aaam_bud nlev")
65
66	if (iim + 1 > 801 .or. jjm + 1 > 401) then
67	print *, ' Problème de dimension dans aaam_bud'
68	stop 1
69	endif
70
71	dlat = pi / jjm
72	dlon = 2 * pi / real(iim)
73
74	oaam = 0.
75	raam = 0.
76	tmou = 0.
77	tsso = 0.
78	tbls = 0.
79
80	! Mountain height, pressure and barotropic wind:
81
82	! North pole values (j = 1):
83
84	ub(1, 1) = 0.
85	vb(1, 1) = 0.
86	do k = 1, nlev
87	ub(1, 1) = ub(1, 1) + u(1, k) * (p(1, k) - p(1, k + 1)) / rg
88	vb(1, 1) = vb(1, 1) + v(1, k) * (p(1, k) - p(1, k + 1)) / rg
89	enddo
90
91	zlat(1) = plat(1) * pi / 180.
92
93	do i = 1, iim + 1
94	zs(i, 1) = phis(1) / rg
95	ps(i, 1) = p(1, 1)
96	ub(i, 1) = ub(1, 1)
97	vb(i, 1) = vb(1, 1)
98	ssou(i, 1) = dragu(1) + liftu(1)
99	ssov(i, 1) = dragv(1) + liftv(1)
100	blsu(i, 1) = phyu(1) - dragu(1) - liftu(1)
101	blsv(i, 1) = phyv(1) - dragv(1) - liftv(1)
102	enddo
103
104	l = 1
105	do j = 2, jjm
106	! Values at Greenwich (Periodicity)
107
108	zs(iim + 1, j) = phis(l + 1) / rg
109	ps(iim + 1, j) = p(l + 1, 1)
110	ssou(iim + 1, j) = dragu(l + 1) + liftu(l + 1)
111	ssov(iim + 1, j) = dragv(l + 1) + liftv(l + 1)
112	blsu(iim + 1, j) = phyu(l + 1) - dragu(l + 1) - liftu(l + 1)
113	blsv(iim + 1, j) = phyv(l + 1) - dragv(l + 1) - liftv(l + 1)
114	zlon(iim + 1) = - plon(l + 1) * pi / 180.
115	zlat(j) = plat(l + 1) * pi / 180.
116
117	ub(iim + 1, j) = 0.
118	vb(iim + 1, j) = 0.
119	do k = 1, nlev
120	ub(iim + 1, j) = ub(iim + 1, j) &
121	+ u(l + 1, k) * (p(l + 1, k) - p(l + 1, k + 1)) / rg
122	vb(iim + 1, j) = vb(iim + 1, j) &
123	+ v(l + 1, k) * (p(l + 1, k) - p(l + 1, k + 1)) / rg
124	enddo
125
126	do i = 1, iim
127	l = l + 1
128	zs(i, j) = phis(l) / rg
129	ps(i, j) = p(l, 1)
130	ssou(i, j) = dragu(l) + liftu(l)
131	ssov(i, j) = dragv(l) + liftv(l)
132	blsu(i, j) = phyu(l) - dragu(l) - liftu(l)
133	blsv(i, j) = phyv(l) - dragv(l) - liftv(l)
134	zlon(i) = plon(l) * pi / 180.
135
136	ub(i, j) = 0.
137	vb(i, j) = 0.
138	do k = 1, nlev
139	ub(i, j) = ub(i, j) + u(l, k) * (p(l, k) - p(l, k + 1)) / rg
140	vb(i, j) = vb(i, j) + v(l, k) * (p(l, k) - p(l, k + 1)) / rg
141	enddo
142	enddo
143	enddo
144
145	! South Pole
146
147	l = l + 1
148	ub(1, jjm + 1) = 0.
149	vb(1, jjm + 1) = 0.
150	do k = 1, nlev
151	ub(1, jjm + 1) = ub(1, jjm + 1) + u(l, k) * (p(l, k) - p(l, k + 1)) / rg
152	vb(1, jjm + 1) = vb(1, jjm + 1) + v(l, k) * (p(l, k) - p(l, k + 1)) / rg
153	enddo
154	zlat(jjm + 1) = plat(l) * pi / 180.
155
156	do i = 1, iim + 1
157	zs(i, jjm + 1) = phis(l) / rg
158	ps(i, jjm + 1) = p(l, 1)
159	ssou(i, jjm + 1) = dragu(l) + liftu(l)
160	ssov(i, jjm + 1) = dragv(l) + liftv(l)
161	blsu(i, jjm + 1) = phyu(l) - dragu(l) - liftu(l)
162	blsv(i, jjm + 1) = phyv(l) - dragv(l) - liftv(l)
163	ub(i, jjm + 1) = ub(1, jjm + 1)
164	vb(i, jjm + 1) = vb(1, jjm + 1)
165	enddo
166
167	! Moment angulaire
168
169	DO j = 1, jjm
170	DO i = 1, iim
171	raam(1) = raam(1) - ra*3 dlon * dlat * 0.5 * (cos(zlon(i )) &
172	* sin(zlat(j )) * cos(zlat(j )) * ub(i , j ) + cos(zlon(i )) &
173	* sin(zlat(j + 1)) * cos(zlat(j + 1)) * ub(i , j + 1)) &
174	+ ra*3 dlon * dlat * 0.5 * (sin(zlon(i )) * cos(zlat(j )) &
175	* vb(i , j ) + sin(zlon(i )) * cos(zlat(j + 1)) * vb(i , j + 1))
176
177	oaam(1) = oaam(1) - ome * ra*4 dlon * dlat / rg * 0.5 &
178	* (cos(zlon(i )) * cos(zlat(j ))*2 sin(zlat(j )) &
179	* ps(i , j ) + cos(zlon(i )) * cos(zlat(j + 1))**2 &
180	* sin(zlat(j + 1)) * ps(i , j + 1))
181
182	raam(2) = raam(2) - ra*3 dlon * dlat * 0.5 * (sin(zlon(i )) &
183	* sin(zlat(j )) * cos(zlat(j )) * ub(i , j ) + sin(zlon(i )) &
184	* sin(zlat(j + 1)) * cos(zlat(j + 1)) * ub(i , j + 1)) &
185	- ra*3 dlon * dlat * 0.5 * (cos(zlon(i )) * cos(zlat(j )) &
186	* vb(i , j ) + cos(zlon(i )) * cos(zlat(j + 1)) * vb(i , j + 1))
187
188	oaam(2) = oaam(2) - ome * ra*4 dlon * dlat / rg * 0.5 &
189	* (sin(zlon(i )) * cos(zlat(j ))*2 sin(zlat(j )) &
190	* ps(i , j ) + sin(zlon(i )) * cos(zlat(j + 1))**2 &
191	* sin(zlat(j + 1)) * ps(i , j + 1))
192
193	raam(3) = raam(3) + ra*3 dlon * dlat * 0.5 * (cos(zlat(j))**2 &
194	* ub(i, j) + cos(zlat(j + 1))*2 ub(i, j + 1))
195
196	oaam(3) = oaam(3) + ome * ra*4 dlon * dlat / rg * 0.5 &
197	* (cos(zlat(j))*3 ps(i, j) + cos(zlat(j + 1))**3 &
198	* ps(i, j + 1))
199	ENDDO
200	ENDDO
201
202	! Couple des montagnes :
203
204	DO j = 1, jjm
205	DO i = 1, iim
206	tmou(1) = tmou(1) - ra*2 dlon * 0.5 * sin(zlon(i)) &
207	* (zs(i, j) - zs(i, j + 1)) &
208	* (cos(zlat(j + 1)) * ps(i, j + 1) + cos(zlat(j)) * ps(i, j))
209	tmou(2) = tmou(2) + ra*2 dlon * 0.5 * cos(zlon(i)) &
210	* (zs(i, j) - zs(i, j + 1)) &
211	* (cos(zlat(j + 1)) * ps(i, j + 1) + cos(zlat(j)) * ps(i, j))
212	ENDDO
213	ENDDO
214
215	DO j = 2, jjm
216	DO i = 1, iim
217	tmou(1) = tmou(1) + ra*2 dlat * 0.5 * sin(zlat(j)) &
218	* (zs(i + 1, j) - zs(i, j)) &
219	* (cos(zlon(i + 1)) * ps(i + 1, j) + cos(zlon(i)) * ps(i, j))
220	tmou(2) = tmou(2) + ra*2 dlat * 0.5 * sin(zlat(j)) &
221	* (zs(i + 1, j) - zs(i, j)) &
222	* (sin(zlon(i + 1)) * ps(i + 1, j) + sin(zlon(i)) * ps(i, j))
223	tmou(3) = tmou(3) - ra*2 dlat * 0.5* cos(zlat(j)) &
224	* (zs(i + 1, j) - zs(i, j)) * (ps(i + 1, j) + ps(i, j))
225	ENDDO
226	ENDDO
227
228	! Couples des differentes friction au sol :
229
230	DO j = 2, jjm
231	DO i = 1, iim
232	tsso(1) = tsso(1) - ra*3 cos(zlat(j)) * dlon * dlat* &
233	ssou(i, j) * sin(zlat(j)) * cos(zlon(i)) &
234	+ ra*3 cos(zlat(j)) * dlon * dlat* &
235	ssov(i, j) * sin(zlon(i))
236
237	tsso(2) = tsso(2) - ra*3 cos(zlat(j)) * dlon * dlat* &
238	ssou(i, j) * sin(zlat(j)) * sin(zlon(i)) &
239	- ra*3 cos(zlat(j)) * dlon * dlat* &
240	ssov(i, j) * cos(zlon(i))
241
242	tsso(3) = tsso(3) + ra*3 cos(zlat(j)) * dlon * dlat* &
243	ssou(i, j) * cos(zlat(j))
244
245	tbls(1) = tbls(1) - ra*3 cos(zlat(j)) * dlon * dlat* &
246	blsu(i, j) * sin(zlat(j)) * cos(zlon(i)) &
247	+ ra*3 cos(zlat(j)) * dlon * dlat* &
248	blsv(i, j) * sin(zlon(i))
249
250	tbls(2) = tbls(2) - ra*3 cos(zlat(j)) * dlon * dlat* &
251	blsu(i, j) * sin(zlat(j)) * sin(zlon(i)) &
252	- ra*3 cos(zlat(j)) * dlon * dlat* &
253	blsv(i, j) * cos(zlon(i))
254
255	tbls(3) = tbls(3) + ra*3 cos(zlat(j)) * dlon * dlat* &
256	blsu(i, j) * cos(zlat(j))
257	ENDDO
258	ENDDO
259
260	aam = raam(3)
261	torsfc = tmou(3) + tsso(3) + tbls(3)
262
263	END subroutine aaam_bud
264
265	end module aaam_bud_m