Sources/dyn3d/pentes_ini.f


! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/pentes_ini.F,v 1.1.1.1 2004/05/19
! 12:53:07 lmdzadmin Exp $

SUBROUTINE pentes_ini(q, w, masse, pbaru, pbarv, mode)

  USE comconst
  USE dynetat0_m, only: rlonv, rlonu
  USE dimens_m
  USE disvert_m
  USE nr_util, ONLY: pi
  USE paramet_m

  IMPLICIT NONE

  ! =======================================================================
  ! Adaptation LMDZ:  A.Armengaud (LGGE)
  ! ----------------

  ! ********************************************************************
  ! Transport des traceurs par la methode des pentes
  ! ********************************************************************
  ! Reference possible : Russel. G.L., Lerner J.A.:
  ! A new Finite-Differencing Scheme for Traceur Transport
  ! Equation , Journal of Applied Meteorology, pp 1483-1498,dec. 81
  ! ********************************************************************
  ! q,w,masse,pbaru et pbarv
  ! sont des arguments d'entree  pour le s-pg ....

  ! =======================================================================


  ! Arguments:
  ! ----------
  INTEGER mode
  REAL, INTENT (IN) :: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
  REAL q(iip1, jjp1, llm, 0:3)
  REAL w(ip1jmp1, llm)
  REAL masse(iip1, jjp1, llm)
  ! Local:
  ! ------
  LOGICAL limit
  REAL sm(iip1, jjp1, llm)
  REAL s0(iip1, jjp1, llm), sx(iip1, jjp1, llm)
  REAL sy(iip1, jjp1, llm), sz(iip1, jjp1, llm)
  REAL masn, mass, zz
  INTEGER i, j, l, iq

  ! modif Fred 24 03 96

  REAL sinlon(iip1), sinlondlon(iip1)
  REAL coslon(iip1), coslondlon(iip1)
  SAVE sinlon, coslon, sinlondlon, coslondlon
  REAL dyn1, dyn2, dys1, dys2
  REAL qpn, qps, dqzpn, dqzps
  REAL smn, sms, s0n, s0s, sxn(iip1), sxs(iip1)
  REAL qmin, pente_max

  REAL ssum
  INTEGER ismax, ismin, lati, latf
  EXTERNAL ssum, convflu, ismin, ismax
  LOGICAL first
  SAVE first
  ! fin modif

  ! EXTERNAL masskg
  EXTERNAL advx
  EXTERNAL advy
  EXTERNAL advz

  ! modif Fred 24 03 96
  DATA first/.TRUE./

  limit = .TRUE.
  pente_max = 2
  ! if (mode.eq.1.or.mode.eq.3) then
  ! if (mode.eq.1) then
  IF (mode>=1) THEN
    lati = 2
    latf = jjm
  ELSE
    lati = 1
    latf = jjp1
  END IF

  qmin = 0.4995
  qmin = 0.
  IF (first) THEN
    PRINT *, 'SCHEMA AMONT NOUVEAU'
    first = .FALSE.
    DO i = 2, iip1
      coslon(i) = cos(rlonv(i))
      sinlon(i) = sin(rlonv(i))
      coslondlon(i) = coslon(i)*(rlonu(i)-rlonu(i-1))/pi
      sinlondlon(i) = sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
      PRINT *, coslondlon(i), sinlondlon(i)
    END DO
    coslon(1) = coslon(iip1)
    coslondlon(1) = coslondlon(iip1)
    sinlon(1) = sinlon(iip1)
    sinlondlon(1) = sinlondlon(iip1)
    PRINT *, 'sum sinlondlon ', ssum(iim, sinlondlon, 1)/sinlondlon(1)
    PRINT *, 'sum coslondlon ', ssum(iim, coslondlon, 1)/coslondlon(1)
    DO l = 1, llm
      DO j = 1, jjp1
        DO i = 1, iip1
          q(i, j, l, 1) = 0.
          q(i, j, l, 2) = 0.
          q(i, j, l, 3) = 0.
        END DO
      END DO
    END DO

  END IF

  ! *** q contient les qqtes de traceur avant l'advection

  ! *** Affectation des tableaux S a partir de Q
  ! *** Rem : utilisation de SCOPY ulterieurement

  DO l = 1, llm
    DO j = 1, jjp1
      DO i = 1, iip1
        s0(i, j, llm+1-l) = q(i, j, l, 0)
        sx(i, j, llm+1-l) = q(i, j, l, 1)
        sy(i, j, llm+1-l) = q(i, j, l, 2)
        sz(i, j, llm+1-l) = q(i, j, l, 3)
      END DO
    END DO
  END DO

  ! *** On calcule la masse d'air en kg

  DO l = 1, llm
    DO j = 1, jjp1
      DO i = 1, iip1
        sm(i, j, llm+1-l) = masse(i, j, l)
      END DO
    END DO
  END DO

  ! *** On converti les champs S en atome (resp. kg)
  ! *** Les routines d'advection traitent les champs
  ! *** a advecter si ces derniers sont en atome (resp. kg)
  ! *** A optimiser !!!

  DO l = 1, llm
    DO j = 1, jjp1
      DO i = 1, iip1
        s0(i, j, l) = s0(i, j, l)*sm(i, j, l)
        sx(i, j, l) = sx(i, j, l)*sm(i, j, l)
        sy(i, j, l) = sy(i, j, l)*sm(i, j, l)
        sz(i, j, l) = sz(i, j, l)*sm(i, j, l)
      END DO
    END DO
  END DO

  ! *** Appel des subroutines d'advection en X, en Y et en Z
  ! *** Advection avec "time-splitting"

  IF (mode==2) THEN
    DO l = 1, llm
      s0s = 0.
      s0n = 0.
      dyn1 = 0.
      dys1 = 0.
      dyn2 = 0.
      dys2 = 0.
      smn = 0.
      sms = 0.
      DO i = 1, iim
        smn = smn + sm(i, 1, l)
        sms = sms + sm(i, jjp1, l)
        s0n = s0n + s0(i, 1, l)
        s0s = s0s + s0(i, jjp1, l)
        zz = sy(i, 1, l)/sm(i, 1, l)
        dyn1 = dyn1 + sinlondlon(i)*zz
        dyn2 = dyn2 + coslondlon(i)*zz
        zz = sy(i, jjp1, l)/sm(i, jjp1, l)
        dys1 = dys1 + sinlondlon(i)*zz
        dys2 = dys2 + coslondlon(i)*zz
      END DO
      DO i = 1, iim
        sy(i, 1, l) = dyn1*sinlon(i) + dyn2*coslon(i)
        sy(i, jjp1, l) = dys1*sinlon(i) + dys2*coslon(i)
      END DO
      DO i = 1, iim
        s0(i, 1, l) = s0n/smn + sy(i, 1, l)
        s0(i, jjp1, l) = s0s/sms - sy(i, jjp1, l)
      END DO

      s0(iip1, 1, l) = s0(1, 1, l)
      s0(iip1, jjp1, l) = s0(1, jjp1, l)

      DO i = 1, iim
        sxn(i) = s0(i+1, 1, l) - s0(i, 1, l)
        sxs(i) = s0(i+1, jjp1, l) - s0(i, jjp1, l)
        ! on rerentre les masses
      END DO
      DO i = 1, iim
        sy(i, 1, l) = sy(i, 1, l)*sm(i, 1, l)
        sy(i, jjp1, l) = sy(i, jjp1, l)*sm(i, jjp1, l)
        s0(i, 1, l) = s0(i, 1, l)*sm(i, 1, l)
        s0(i, jjp1, l) = s0(i, jjp1, l)*sm(i, jjp1, l)
      END DO
      sxn(iip1) = sxn(1)
      sxs(iip1) = sxs(1)
      DO i = 1, iim
        sx(i+1, 1, l) = 0.25*(sxn(i)+sxn(i+1))*sm(i+1, 1, l)
        sx(i+1, jjp1, l) = 0.25*(sxs(i)+sxs(i+1))*sm(i+1, jjp1, l)
      END DO
      s0(iip1, 1, l) = s0(1, 1, l)
      s0(iip1, jjp1, l) = s0(1, jjp1, l)
      sy(iip1, 1, l) = sy(1, 1, l)
      sy(iip1, jjp1, l) = sy(1, jjp1, l)
      sx(1, 1, l) = sx(iip1, 1, l)
      sx(1, jjp1, l) = sx(iip1, jjp1, l)
    END DO
  END IF

  IF (mode==4) THEN
    DO l = 1, llm
      DO i = 1, iip1
        sx(i, 1, l) = 0.
        sx(i, jjp1, l) = 0.
        sy(i, 1, l) = 0.
        sy(i, jjp1, l) = 0.
      END DO
    END DO
  END IF
  CALL limx(s0, sx, sm, pente_max)
  CALL advx(limit, .5*dtvr, pbaru, sm, s0, sx, sy, sz, lati, latf)
  IF (mode==4) THEN
    DO l = 1, llm
      DO i = 1, iip1
        sx(i, 1, l) = 0.
        sx(i, jjp1, l) = 0.
        sy(i, 1, l) = 0.
        sy(i, jjp1, l) = 0.
      END DO
    END DO
  END IF
  CALL limy(s0, sy, sm, pente_max)
  CALL advy(limit, .5*dtvr, pbarv, sm, s0, sx, sy, sz)
  DO j = 1, jjp1
    DO i = 1, iip1
      sz(i, j, 1) = 0.
      sz(i, j, llm) = 0.
    END DO
  END DO
  CALL limz(s0, sz, sm, pente_max)
  CALL advz(limit, dtvr, w, sm, s0, sx, sy, sz)
  IF (mode==4) THEN
    DO l = 1, llm
      DO i = 1, iip1
        sx(i, 1, l) = 0.
        sx(i, jjp1, l) = 0.
        sy(i, 1, l) = 0.
        sy(i, jjp1, l) = 0.
      END DO
    END DO
  END IF
  CALL limy(s0, sy, sm, pente_max)
  CALL advy(limit, .5*dtvr, pbarv, sm, s0, sx, sy, sz)
  DO l = 1, llm
    DO j = 1, jjp1
      sm(iip1, j, l) = sm(1, j, l)
      s0(iip1, j, l) = s0(1, j, l)
      sx(iip1, j, l) = sx(1, j, l)
      sy(iip1, j, l) = sy(1, j, l)
      sz(iip1, j, l) = sz(1, j, l)
    END DO
  END DO


  IF (mode==4) THEN
    DO l = 1, llm
      DO i = 1, iip1
        sx(i, 1, l) = 0.
        sx(i, jjp1, l) = 0.
        sy(i, 1, l) = 0.
        sy(i, jjp1, l) = 0.
      END DO
    END DO
  END IF
  CALL limx(s0, sx, sm, pente_max)
  CALL advx(limit, .5*dtvr, pbaru, sm, s0, sx, sy, sz, lati, latf)
  ! ***   On repasse les S dans la variable q directement 14/10/94
  ! On revient a des rapports de melange en divisant par la masse

  ! En dehors des poles:

  DO l = 1, llm
    DO j = 1, jjp1
      DO i = 1, iim
        q(i, j, llm+1-l, 0) = s0(i, j, l)/sm(i, j, l)
        q(i, j, llm+1-l, 1) = sx(i, j, l)/sm(i, j, l)
        q(i, j, llm+1-l, 2) = sy(i, j, l)/sm(i, j, l)
        q(i, j, llm+1-l, 3) = sz(i, j, l)/sm(i, j, l)
      END DO
    END DO
  END DO

  ! Traitements specifiques au pole

  IF (mode>=1) THEN
    DO l = 1, llm
      ! filtrages aux poles
      masn = ssum(iim, sm(1,1,l), 1)
      mass = ssum(iim, sm(1,jjp1,l), 1)
      qpn = ssum(iim, s0(1,1,l), 1)/masn
      qps = ssum(iim, s0(1,jjp1,l), 1)/mass
      dqzpn = ssum(iim, sz(1,1,l), 1)/masn
      dqzps = ssum(iim, sz(1,jjp1,l), 1)/mass
      DO i = 1, iip1
        q(i, 1, llm+1-l, 3) = dqzpn
        q(i, jjp1, llm+1-l, 3) = dqzps
        q(i, 1, llm+1-l, 0) = qpn
        q(i, jjp1, llm+1-l, 0) = qps
      END DO
      IF (mode==3) THEN
        dyn1 = 0.
        dys1 = 0.
        dyn2 = 0.
        dys2 = 0.
        DO i = 1, iim
          dyn1 = dyn1 + sinlondlon(i)*sy(i, 1, l)/sm(i, 1, l)
          dyn2 = dyn2 + coslondlon(i)*sy(i, 1, l)/sm(i, 1, l)
          dys1 = dys1 + sinlondlon(i)*sy(i, jjp1, l)/sm(i, jjp1, l)
          dys2 = dys2 + coslondlon(i)*sy(i, jjp1, l)/sm(i, jjp1, l)
        END DO
        DO i = 1, iim
          q(i, 1, llm+1-l, 2) = (sinlon(i)*dyn1+coslon(i)*dyn2)
          q(i, 1, llm+1-l, 0) = q(i, 1, llm+1-l, 0) + q(i, 1, llm+1-l, 2)
          q(i, jjp1, llm+1-l, 2) = (sinlon(i)*dys1+coslon(i)*dys2)
          q(i, jjp1, llm+1-l, 0) = q(i, jjp1, llm+1-l, 0) - &
            q(i, jjp1, llm+1-l, 2)
        END DO
      END IF
      IF (mode==1) THEN
        ! on filtre les valeurs au bord de la "grande maille pole"
        dyn1 = 0.
        dys1 = 0.
        dyn2 = 0.
        dys2 = 0.
        DO i = 1, iim
          zz = s0(i, 2, l)/sm(i, 2, l) - q(i, 1, llm+1-l, 0)
          dyn1 = dyn1 + sinlondlon(i)*zz
          dyn2 = dyn2 + coslondlon(i)*zz
          zz = q(i, jjp1, llm+1-l, 0) - s0(i, jjm, l)/sm(i, jjm, l)
          dys1 = dys1 + sinlondlon(i)*zz
          dys2 = dys2 + coslondlon(i)*zz
        END DO
        DO i = 1, iim
          q(i, 1, llm+1-l, 2) = (sinlon(i)*dyn1+coslon(i)*dyn2)/2.
          q(i, 1, llm+1-l, 0) = q(i, 1, llm+1-l, 0) + q(i, 1, llm+1-l, 2)
          q(i, jjp1, llm+1-l, 2) = (sinlon(i)*dys1+coslon(i)*dys2)/2.
          q(i, jjp1, llm+1-l, 0) = q(i, jjp1, llm+1-l, 0) - &
            q(i, jjp1, llm+1-l, 2)
        END DO
        q(iip1, 1, llm+1-l, 0) = q(1, 1, llm+1-l, 0)
        q(iip1, jjp1, llm+1-l, 0) = q(1, jjp1, llm+1-l, 0)

        DO i = 1, iim
          sxn(i) = q(i+1, 1, llm+1-l, 0) - q(i, 1, llm+1-l, 0)
          sxs(i) = q(i+1, jjp1, llm+1-l, 0) - q(i, jjp1, llm+1-l, 0)
        END DO
        sxn(iip1) = sxn(1)
        sxs(iip1) = sxs(1)
        DO i = 1, iim
          q(i+1, 1, llm+1-l, 1) = 0.25*(sxn(i)+sxn(i+1))
          q(i+1, jjp1, llm+1-l, 1) = 0.25*(sxs(i)+sxs(i+1))
        END DO
        q(1, 1, llm+1-l, 1) = q(iip1, 1, llm+1-l, 1)
        q(1, jjp1, llm+1-l, 1) = q(iip1, jjp1, llm+1-l, 1)

      END IF

    END DO
  END IF

  ! bouclage en longitude
  DO iq = 0, 3
    DO l = 1, llm
      DO j = 1, jjp1
        q(iip1, j, l, iq) = q(1, j, l, iq)
      END DO
    END DO
  END DO

  DO l = 1, llm
    DO j = 1, jjp1
      DO i = 1, iip1
        IF (q(i,j,l,0)<0.) THEN
          q(i, j, l, 0) = 0.
        END IF
      END DO
    END DO
  END DO

  DO l = 1, llm
    DO j = 1, jjp1
      DO i = 1, iip1
        IF (q(i,j,l,0)<qmin) PRINT *, 'apres pentes, s0(', i, ',', j, ',', l, &
          ')=', q(i, j, l, 0)
      END DO
    END DO
  END DO
  RETURN
END SUBROUTINE pentes_ini
1	guez	3
2	guez	81	! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/pentes_ini.F,v 1.1.1.1 2004/05/19
3			! 12:53:07 lmdzadmin Exp $
4	guez	3
5	guez	81	SUBROUTINE pentes_ini(q, w, masse, pbaru, pbarv, mode)
6	guez	139
7			USE comconst
8			USE dynetat0_m, only: rlonv, rlonu
9	guez	81	USE dimens_m
10			USE disvert_m
11			USE nr_util, ONLY: pi
12	guez	139	USE paramet_m
13
14	guez	81	IMPLICIT NONE
15	guez	3
16	guez	81	! =======================================================================
17			! Adaptation LMDZ: A.Armengaud (LGGE)
18			! ----------------
19	guez	3
20	guez	81	! ********************************************************************
21			! Transport des traceurs par la methode des pentes
22			! ********************************************************************
23			! Reference possible : Russel. G.L., Lerner J.A.:
24			! A new Finite-Differencing Scheme for Traceur Transport
25			! Equation , Journal of Applied Meteorology, pp 1483-1498,dec. 81
26			! ********************************************************************
27			! q,w,masse,pbaru et pbarv
28			! sont des arguments d'entree pour le s-pg ....
29	guez	3
30	guez	81	! =======================================================================
31	guez	3
32
33
34	guez	81	! Arguments:
35			! ----------
36			INTEGER mode
37			REAL, INTENT (IN) :: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
38			REAL q(iip1, jjp1, llm, 0:3)
39			REAL w(ip1jmp1, llm)
40			REAL masse(iip1, jjp1, llm)
41			! Local:
42			! ------
43			LOGICAL limit
44			REAL sm(iip1, jjp1, llm)
45			REAL s0(iip1, jjp1, llm), sx(iip1, jjp1, llm)
46			REAL sy(iip1, jjp1, llm), sz(iip1, jjp1, llm)
47			REAL masn, mass, zz
48			INTEGER i, j, l, iq
49	guez	3
50	guez	81	! modif Fred 24 03 96
51	guez	3
52	guez	81	REAL sinlon(iip1), sinlondlon(iip1)
53			REAL coslon(iip1), coslondlon(iip1)
54			SAVE sinlon, coslon, sinlondlon, coslondlon
55			REAL dyn1, dyn2, dys1, dys2
56			REAL qpn, qps, dqzpn, dqzps
57			REAL smn, sms, s0n, s0s, sxn(iip1), sxs(iip1)
58	guez	105	REAL qmin, pente_max
59	guez	3
60	guez	81	REAL ssum
61			INTEGER ismax, ismin, lati, latf
62			EXTERNAL ssum, convflu, ismin, ismax
63			LOGICAL first
64			SAVE first
65			! fin modif
66	guez	3
67	guez	81	! EXTERNAL masskg
68			EXTERNAL advx
69			EXTERNAL advy
70			EXTERNAL advz
71	guez	3
72	guez	81	! modif Fred 24 03 96
73			DATA first/.TRUE./
74	guez	3
75	guez	81	limit = .TRUE.
76			pente_max = 2
77			! if (mode.eq.1.or.mode.eq.3) then
78			! if (mode.eq.1) then
79			IF (mode>=1) THEN
80			lati = 2
81			latf = jjm
82			ELSE
83			lati = 1
84			latf = jjp1
85			END IF
86	guez	3
87	guez	81	qmin = 0.4995
88			qmin = 0.
89			IF (first) THEN
90			PRINT *, 'SCHEMA AMONT NOUVEAU'
91			first = .FALSE.
92			DO i = 2, iip1
93			coslon(i) = cos(rlonv(i))
94			sinlon(i) = sin(rlonv(i))
95			coslondlon(i) = coslon(i)*(rlonu(i)-rlonu(i-1))/pi
96			sinlondlon(i) = sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
97			PRINT *, coslondlon(i), sinlondlon(i)
98			END DO
99			coslon(1) = coslon(iip1)
100			coslondlon(1) = coslondlon(iip1)
101			sinlon(1) = sinlon(iip1)
102			sinlondlon(1) = sinlondlon(iip1)
103			PRINT *, 'sum sinlondlon ', ssum(iim, sinlondlon, 1)/sinlondlon(1)
104			PRINT *, 'sum coslondlon ', ssum(iim, coslondlon, 1)/coslondlon(1)
105			DO l = 1, llm
106			DO j = 1, jjp1
107			DO i = 1, iip1
108			q(i, j, l, 1) = 0.
109			q(i, j, l, 2) = 0.
110			q(i, j, l, 3) = 0.
111			END DO
112			END DO
113			END DO
114	guez	3
115	guez	81	END IF
116	guez	3
117	guez	81	! *** q contient les qqtes de traceur avant l'advection
118	guez	3
119	guez	81	! *** Affectation des tableaux S a partir de Q
120			! *** Rem : utilisation de SCOPY ulterieurement
121	guez	3
122	guez	81	DO l = 1, llm
123			DO j = 1, jjp1
124			DO i = 1, iip1
125			s0(i, j, llm+1-l) = q(i, j, l, 0)
126			sx(i, j, llm+1-l) = q(i, j, l, 1)
127			sy(i, j, llm+1-l) = q(i, j, l, 2)
128			sz(i, j, llm+1-l) = q(i, j, l, 3)
129			END DO
130			END DO
131			END DO
132	guez	3
133	guez	81	! *** On calcule la masse d'air en kg
134	guez	3
135	guez	81	DO l = 1, llm
136			DO j = 1, jjp1
137			DO i = 1, iip1
138			sm(i, j, llm+1-l) = masse(i, j, l)
139			END DO
140			END DO
141			END DO
142	guez	3
143	guez	81	! *** On converti les champs S en atome (resp. kg)
144			! *** Les routines d'advection traitent les champs
145			! *** a advecter si ces derniers sont en atome (resp. kg)
146			! *** A optimiser !!!
147	guez	3
148	guez	81	DO l = 1, llm
149			DO j = 1, jjp1
150			DO i = 1, iip1
151			s0(i, j, l) = s0(i, j, l)*sm(i, j, l)
152			sx(i, j, l) = sx(i, j, l)*sm(i, j, l)
153			sy(i, j, l) = sy(i, j, l)*sm(i, j, l)
154			sz(i, j, l) = sz(i, j, l)*sm(i, j, l)
155			END DO
156			END DO
157			END DO
158	guez	3
159	guez	81	! *** Appel des subroutines d'advection en X, en Y et en Z
160			! *** Advection avec "time-splitting"
161	guez	3
162	guez	81	IF (mode==2) THEN
163			DO l = 1, llm
164			s0s = 0.
165			s0n = 0.
166			dyn1 = 0.
167			dys1 = 0.
168			dyn2 = 0.
169			dys2 = 0.
170			smn = 0.
171			sms = 0.
172			DO i = 1, iim
173			smn = smn + sm(i, 1, l)
174			sms = sms + sm(i, jjp1, l)
175			s0n = s0n + s0(i, 1, l)
176			s0s = s0s + s0(i, jjp1, l)
177			zz = sy(i, 1, l)/sm(i, 1, l)
178			dyn1 = dyn1 + sinlondlon(i)*zz
179			dyn2 = dyn2 + coslondlon(i)*zz
180			zz = sy(i, jjp1, l)/sm(i, jjp1, l)
181			dys1 = dys1 + sinlondlon(i)*zz
182			dys2 = dys2 + coslondlon(i)*zz
183			END DO
184			DO i = 1, iim
185			sy(i, 1, l) = dyn1sinlon(i) + dyn2coslon(i)
186			sy(i, jjp1, l) = dys1sinlon(i) + dys2coslon(i)
187			END DO
188			DO i = 1, iim
189			s0(i, 1, l) = s0n/smn + sy(i, 1, l)
190			s0(i, jjp1, l) = s0s/sms - sy(i, jjp1, l)
191			END DO
192	guez	3
193	guez	81	s0(iip1, 1, l) = s0(1, 1, l)
194			s0(iip1, jjp1, l) = s0(1, jjp1, l)
195	guez	3
196	guez	81	DO i = 1, iim
197			sxn(i) = s0(i+1, 1, l) - s0(i, 1, l)
198			sxs(i) = s0(i+1, jjp1, l) - s0(i, jjp1, l)
199			! on rerentre les masses
200			END DO
201			DO i = 1, iim
202			sy(i, 1, l) = sy(i, 1, l)*sm(i, 1, l)
203			sy(i, jjp1, l) = sy(i, jjp1, l)*sm(i, jjp1, l)
204			s0(i, 1, l) = s0(i, 1, l)*sm(i, 1, l)
205			s0(i, jjp1, l) = s0(i, jjp1, l)*sm(i, jjp1, l)
206			END DO
207			sxn(iip1) = sxn(1)
208			sxs(iip1) = sxs(1)
209			DO i = 1, iim
210			sx(i+1, 1, l) = 0.25(sxn(i)+sxn(i+1))sm(i+1, 1, l)
211			sx(i+1, jjp1, l) = 0.25(sxs(i)+sxs(i+1))sm(i+1, jjp1, l)
212			END DO
213			s0(iip1, 1, l) = s0(1, 1, l)
214			s0(iip1, jjp1, l) = s0(1, jjp1, l)
215			sy(iip1, 1, l) = sy(1, 1, l)
216			sy(iip1, jjp1, l) = sy(1, jjp1, l)
217			sx(1, 1, l) = sx(iip1, 1, l)
218			sx(1, jjp1, l) = sx(iip1, jjp1, l)
219			END DO
220			END IF
221	guez	3
222	guez	81	IF (mode==4) THEN
223			DO l = 1, llm
224			DO i = 1, iip1
225			sx(i, 1, l) = 0.
226			sx(i, jjp1, l) = 0.
227			sy(i, 1, l) = 0.
228			sy(i, jjp1, l) = 0.
229			END DO
230			END DO
231			END IF
232			CALL limx(s0, sx, sm, pente_max)
233			CALL advx(limit, .5*dtvr, pbaru, sm, s0, sx, sy, sz, lati, latf)
234			IF (mode==4) THEN
235			DO l = 1, llm
236			DO i = 1, iip1
237			sx(i, 1, l) = 0.
238			sx(i, jjp1, l) = 0.
239			sy(i, 1, l) = 0.
240			sy(i, jjp1, l) = 0.
241			END DO
242			END DO
243			END IF
244			CALL limy(s0, sy, sm, pente_max)
245			CALL advy(limit, .5*dtvr, pbarv, sm, s0, sx, sy, sz)
246			DO j = 1, jjp1
247			DO i = 1, iip1
248			sz(i, j, 1) = 0.
249			sz(i, j, llm) = 0.
250			END DO
251			END DO
252			CALL limz(s0, sz, sm, pente_max)
253			CALL advz(limit, dtvr, w, sm, s0, sx, sy, sz)
254			IF (mode==4) THEN
255			DO l = 1, llm
256			DO i = 1, iip1
257			sx(i, 1, l) = 0.
258			sx(i, jjp1, l) = 0.
259			sy(i, 1, l) = 0.
260			sy(i, jjp1, l) = 0.
261			END DO
262			END DO
263			END IF
264			CALL limy(s0, sy, sm, pente_max)
265			CALL advy(limit, .5*dtvr, pbarv, sm, s0, sx, sy, sz)
266			DO l = 1, llm
267			DO j = 1, jjp1
268			sm(iip1, j, l) = sm(1, j, l)
269			s0(iip1, j, l) = s0(1, j, l)
270			sx(iip1, j, l) = sx(1, j, l)
271			sy(iip1, j, l) = sy(1, j, l)
272			sz(iip1, j, l) = sz(1, j, l)
273			END DO
274			END DO
275	guez	3
276
277	guez	81	IF (mode==4) THEN
278			DO l = 1, llm
279			DO i = 1, iip1
280			sx(i, 1, l) = 0.
281			sx(i, jjp1, l) = 0.
282			sy(i, 1, l) = 0.
283			sy(i, jjp1, l) = 0.
284			END DO
285			END DO
286			END IF
287			CALL limx(s0, sx, sm, pente_max)
288			CALL advx(limit, .5*dtvr, pbaru, sm, s0, sx, sy, sz, lati, latf)
289			! *** On repasse les S dans la variable q directement 14/10/94
290			! On revient a des rapports de melange en divisant par la masse
291	guez	3
292	guez	81	! En dehors des poles:
293	guez	3
294	guez	81	DO l = 1, llm
295			DO j = 1, jjp1
296			DO i = 1, iim
297			q(i, j, llm+1-l, 0) = s0(i, j, l)/sm(i, j, l)
298			q(i, j, llm+1-l, 1) = sx(i, j, l)/sm(i, j, l)
299			q(i, j, llm+1-l, 2) = sy(i, j, l)/sm(i, j, l)
300			q(i, j, llm+1-l, 3) = sz(i, j, l)/sm(i, j, l)
301			END DO
302			END DO
303			END DO
304
305			! Traitements specifiques au pole
306
307			IF (mode>=1) THEN
308			DO l = 1, llm
309			! filtrages aux poles
310			masn = ssum(iim, sm(1,1,l), 1)
311			mass = ssum(iim, sm(1,jjp1,l), 1)
312			qpn = ssum(iim, s0(1,1,l), 1)/masn
313			qps = ssum(iim, s0(1,jjp1,l), 1)/mass
314			dqzpn = ssum(iim, sz(1,1,l), 1)/masn
315			dqzps = ssum(iim, sz(1,jjp1,l), 1)/mass
316			DO i = 1, iip1
317			q(i, 1, llm+1-l, 3) = dqzpn
318			q(i, jjp1, llm+1-l, 3) = dqzps
319			q(i, 1, llm+1-l, 0) = qpn
320			q(i, jjp1, llm+1-l, 0) = qps
321			END DO
322			IF (mode==3) THEN
323			dyn1 = 0.
324			dys1 = 0.
325			dyn2 = 0.
326			dys2 = 0.
327			DO i = 1, iim
328			dyn1 = dyn1 + sinlondlon(i)*sy(i, 1, l)/sm(i, 1, l)
329			dyn2 = dyn2 + coslondlon(i)*sy(i, 1, l)/sm(i, 1, l)
330			dys1 = dys1 + sinlondlon(i)*sy(i, jjp1, l)/sm(i, jjp1, l)
331			dys2 = dys2 + coslondlon(i)*sy(i, jjp1, l)/sm(i, jjp1, l)
332			END DO
333			DO i = 1, iim
334			q(i, 1, llm+1-l, 2) = (sinlon(i)dyn1+coslon(i)dyn2)
335			q(i, 1, llm+1-l, 0) = q(i, 1, llm+1-l, 0) + q(i, 1, llm+1-l, 2)
336			q(i, jjp1, llm+1-l, 2) = (sinlon(i)dys1+coslon(i)dys2)
337			q(i, jjp1, llm+1-l, 0) = q(i, jjp1, llm+1-l, 0) - &
338			q(i, jjp1, llm+1-l, 2)
339			END DO
340			END IF
341			IF (mode==1) THEN
342			! on filtre les valeurs au bord de la "grande maille pole"
343			dyn1 = 0.
344			dys1 = 0.
345			dyn2 = 0.
346			dys2 = 0.
347			DO i = 1, iim
348			zz = s0(i, 2, l)/sm(i, 2, l) - q(i, 1, llm+1-l, 0)
349			dyn1 = dyn1 + sinlondlon(i)*zz
350			dyn2 = dyn2 + coslondlon(i)*zz
351			zz = q(i, jjp1, llm+1-l, 0) - s0(i, jjm, l)/sm(i, jjm, l)
352			dys1 = dys1 + sinlondlon(i)*zz
353			dys2 = dys2 + coslondlon(i)*zz
354			END DO
355			DO i = 1, iim
356			q(i, 1, llm+1-l, 2) = (sinlon(i)dyn1+coslon(i)dyn2)/2.
357			q(i, 1, llm+1-l, 0) = q(i, 1, llm+1-l, 0) + q(i, 1, llm+1-l, 2)
358			q(i, jjp1, llm+1-l, 2) = (sinlon(i)dys1+coslon(i)dys2)/2.
359			q(i, jjp1, llm+1-l, 0) = q(i, jjp1, llm+1-l, 0) - &
360			q(i, jjp1, llm+1-l, 2)
361			END DO
362			q(iip1, 1, llm+1-l, 0) = q(1, 1, llm+1-l, 0)
363			q(iip1, jjp1, llm+1-l, 0) = q(1, jjp1, llm+1-l, 0)
364
365			DO i = 1, iim
366			sxn(i) = q(i+1, 1, llm+1-l, 0) - q(i, 1, llm+1-l, 0)
367			sxs(i) = q(i+1, jjp1, llm+1-l, 0) - q(i, jjp1, llm+1-l, 0)
368			END DO
369			sxn(iip1) = sxn(1)
370			sxs(iip1) = sxs(1)
371			DO i = 1, iim
372			q(i+1, 1, llm+1-l, 1) = 0.25*(sxn(i)+sxn(i+1))
373			q(i+1, jjp1, llm+1-l, 1) = 0.25*(sxs(i)+sxs(i+1))
374			END DO
375			q(1, 1, llm+1-l, 1) = q(iip1, 1, llm+1-l, 1)
376			q(1, jjp1, llm+1-l, 1) = q(iip1, jjp1, llm+1-l, 1)
377
378			END IF
379
380			END DO
381			END IF
382
383			! bouclage en longitude
384			DO iq = 0, 3
385			DO l = 1, llm
386			DO j = 1, jjp1
387			q(iip1, j, l, iq) = q(1, j, l, iq)
388			END DO
389			END DO
390			END DO
391
392			DO l = 1, llm
393			DO j = 1, jjp1
394			DO i = 1, iip1
395			IF (q(i,j,l,0)<0.) THEN
396			q(i, j, l, 0) = 0.
397			END IF
398			END DO
399			END DO
400			END DO
401
402			DO l = 1, llm
403			DO j = 1, jjp1
404			DO i = 1, iip1
405			IF (q(i,j,l,0)<qmin) PRINT *, 'apres pentes, s0(', i, ',', j, ',', l, &
406			')=', q(i, j, l, 0)
407			END DO
408			END DO
409			END DO
410			RETURN
411			END SUBROUTINE pentes_ini