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