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
2	! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/pentes_ini.F,v 1.1.1.1 2004/05/19
3	! 12:53:07 lmdzadmin Exp $
4
5	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
14	! =======================================================================
15	! Adaptation LMDZ: A.Armengaud (LGGE)
16	! ----------------
17
18	! ********************************************************************
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
28	! =======================================================================
29
30
31
32	! 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
48	! modif Fred 24 03 96
49
50	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
58	REAL ssum
59	INTEGER ismax, ismin, lati, latf
60	EXTERNAL ssum, convflu, ismin, ismax
61	LOGICAL first
62	SAVE first
63	! fin modif
64
65	! EXTERNAL masskg
66	EXTERNAL advx
67	EXTERNAL advy
68	EXTERNAL advz
69
70	! modif Fred 24 03 96
71	DATA first/.TRUE./
72
73	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
85	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
113	END IF
114
115	! *** q contient les qqtes de traceur avant l'advection
116
117	! *** Affectation des tableaux S a partir de Q
118	! *** Rem : utilisation de SCOPY ulterieurement
119
120	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
131	! *** On calcule la masse d'air en kg
132
133	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
141	! *** 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
146	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
157	! *** Appel des subroutines d'advection en X, en Y et en Z
158	! *** Advection avec "time-splitting"
159
160	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
191	s0(iip1, 1, l) = s0(1, 1, l)
192	s0(iip1, jjp1, l) = s0(1, jjp1, l)
193
194	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
220	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
274
275	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
290	! En dehors des poles:
291
292	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