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
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
7	USE comconst
8	USE dynetat0_m, only: rlonv, rlonu
9	USE dimens_m
10	USE disvert_m
11	USE nr_util, ONLY: pi
12	USE paramet_m
13
14	IMPLICIT NONE
15
16	! =======================================================================
17	! Adaptation LMDZ: A.Armengaud (LGGE)
18	! ----------------
19
20	! ********************************************************************
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
30	! =======================================================================
31
32
33
34	! 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
50	! modif Fred 24 03 96
51
52	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	REAL qmin, pente_max
59
60	REAL ssum
61	INTEGER ismax, ismin, lati, latf
62	EXTERNAL ssum, convflu, ismin, ismax
63	LOGICAL first
64	SAVE first
65	! fin modif
66
67	! EXTERNAL masskg
68	EXTERNAL advx
69	EXTERNAL advy
70	EXTERNAL advz
71
72	! modif Fred 24 03 96
73	DATA first/.TRUE./
74
75	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
87	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
115	END IF
116
117	! *** q contient les qqtes de traceur avant l'advection
118
119	! *** Affectation des tableaux S a partir de Q
120	! *** Rem : utilisation de SCOPY ulterieurement
121
122	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
133	! *** On calcule la masse d'air en kg
134
135	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
143	! *** 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
148	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
159	! *** Appel des subroutines d'advection en X, en Y et en Z
160	! *** Advection avec "time-splitting"
161
162	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
193	s0(iip1, 1, l) = s0(1, 1, l)
194	s0(iip1, jjp1, l) = s0(1, jjp1, l)
195
196	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
222	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
276
277	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
292	! En dehors des poles:
293
294	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