trunk/phylmd/aaam_bud.f90

module aaam_bud_m

  implicit none

contains

  subroutine aaam_bud(rg, ome, 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 phyetat0_m, only: rlat, rlon
    USE suphec_m, ONLY: ra

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