Sources/dyn3d/prather.f


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

SUBROUTINE prather(q, w, masse, pbaru, pbarv, nt, dt)

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

  IMPLICIT NONE

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

  ! ************************************************
  ! Transport des traceurs par la methode de prather
  ! Ref :

  ! ************************************************
  ! q,w,pext,pbaru et pbarv : arguments d'entree  pour le s-pg

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


  ! Arguments:
  ! ----------
  INTEGER iq, nt
  REAL, INTENT (IN) :: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
  REAL masse(iip1, jjp1, llm)
  REAL q(iip1, jjp1, llm, 0:9)
  REAL w(ip1jmp1, llm)
  INTEGER ordre, ilim

  ! Local:
  ! ------
  LOGICAL limit
  REAL zq(iip1, jjp1, llm)
  REAL sm(iip1, jjp1, llm)
  REAL s0(iip1, jjp1, llm), sx(iip1, jjp1, llm)
  REAL sy(iip1, jjp1, llm), sz(iip1, jjp1, llm)
  REAL sxx(iip1, jjp1, llm)
  REAL sxy(iip1, jjp1, llm)
  REAL sxz(iip1, jjp1, llm)
  REAL syy(iip1, jjp1, llm)
  REAL syz(iip1, jjp1, llm)
  REAL szz(iip1, jjp1, llm), zz
  INTEGER i, j, l, indice
  REAL sxn(iip1), sxs(iip1)

  REAL sinlon(iip1), sinlondlon(iip1)
  REAL coslon(iip1), coslondlon(iip1)
  REAL qmin, qmax
  SAVE qmin, qmax
  SAVE sinlon, coslon, sinlondlon, coslondlon
  REAL dyn1, dyn2, dys1, dys2, qpn, qps, dqzpn, dqzps
  REAL masn, mass

  REAL ssum
  INTEGER ismax, ismin
  EXTERNAL ssum, convflu, ismin, ismax
  LOGICAL first
  SAVE first
  EXTERNAL advxp, advyp, advzp


  DATA first/.TRUE./
  DATA qmin, qmax/ -1.E33, 1.E33/


  ! ==========================================================================
  ! ==========================================================================
  ! MODIFICATION POUR PAS DE TEMPS ADAPTATIF, dtvr remplace par dt
  ! ==========================================================================
  ! ==========================================================================
  REAL dt
  ! ==========================================================================
  limit = .TRUE.

  IF (first) THEN
    PRINT *, 'SCHEMA PRATHER'
    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
    END DO
    coslon(1) = coslon(iip1)
    coslondlon(1) = coslondlon(iip1)
    sinlon(1) = sinlon(iip1)
    sinlondlon(1) = sinlondlon(iip1)

    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.
          q(i, j, l, 4) = 0.
          q(i, j, l, 5) = 0.
          q(i, j, l, 6) = 0.
          q(i, j, l, 7) = 0.
          q(i, j, l, 8) = 0.
          q(i, j, l, 9) = 0.
        END DO
      END DO
    END DO
  END IF
  ! Fin modif Fred

  ! *** 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

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

  ! *** Affectation des tableaux S a partir de Q

  DO l = 1, llm
    DO j = 1, jjp1
      DO i = 1, iip1
        s0(i, j, l) = q(i, j, llm+1-l, 0)*sm(i, j, l)
        sx(i, j, l) = q(i, j, llm+1-l, 1)*sm(i, j, l)
        sy(i, j, l) = q(i, j, llm+1-l, 2)*sm(i, j, l)
        sz(i, j, l) = q(i, j, llm+1-l, 3)*sm(i, j, l)
        sxx(i, j, l) = q(i, j, llm+1-l, 4)*sm(i, j, l)
        sxy(i, j, l) = q(i, j, llm+1-l, 5)*sm(i, j, l)
        sxz(i, j, l) = q(i, j, llm+1-l, 6)*sm(i, j, l)
        syy(i, j, l) = q(i, j, llm+1-l, 7)*sm(i, j, l)
        syz(i, j, l) = q(i, j, llm+1-l, 8)*sm(i, j, l)
        szz(i, j, l) = q(i, j, llm+1-l, 9)*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"

  ! -----------------------------------------------------------
  DO indice = 1, nt
    CALL advxp(limit, 0.5*dt, pbaru, sm, s0, sx, sy, sz, sxx, sxy, sxz, syy, &
      syz, szz, 1)
  END DO
  DO l = 1, llm
    DO i = 1, iip1
      sy(i, 1, l) = 0.
      sy(i, jjp1, l) = 0.
    END DO
  END DO
  ! ---------------------------------------------------------
  CALL advyp(limit, .5*dt*nt, pbarv, sm, s0, sx, sy, sz, sxx, sxy, sxz, syy, &
    syz, szz, 1)
  ! ---------------------------------------------------------

  ! ---------------------------------------------------------
  DO j = 1, jjp1
    DO i = 1, iip1
      sz(i, j, 1) = 0.
      sz(i, j, llm) = 0.
      sxz(i, j, 1) = 0.
      sxz(i, j, llm) = 0.
      syz(i, j, 1) = 0.
      syz(i, j, llm) = 0.
      szz(i, j, 1) = 0.
      szz(i, j, llm) = 0.
    END DO
  END DO
  CALL advzp(limit, dt*nt, w, sm, s0, sx, sy, sz, sxx, sxy, sxz, syy, syz, &
    szz, 1)
  DO l = 1, llm
    DO i = 1, iip1
      sy(i, 1, l) = 0.
      sy(i, jjp1, l) = 0.
    END DO
  END DO

  ! ---------------------------------------------------------

  ! ---------------------------------------------------------
  CALL advyp(limit, .5*dt*nt, pbarv, sm, s0, sx, sy, sz, sxx, sxy, sxz, syy, &
    syz, szz, 1)
  ! ---------------------------------------------------------
  DO l = 1, llm
    DO j = 1, jjp1
      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)
      sxx(iip1, j, l) = sxx(1, j, l)
      sxy(iip1, j, l) = sxy(1, j, l)
      sxz(iip1, j, l) = sxz(1, j, l)
      syy(iip1, j, l) = syy(1, j, l)
      syz(iip1, j, l) = syz(1, j, l)
      szz(iip1, j, l) = szz(1, j, l)
    END DO
  END DO
  DO indice = 1, nt
    CALL advxp(limit, 0.5*dt, pbaru, sm, s0, sx, sy, sz, sxx, sxy, sxz, syy, &
      syz, szz, 1)
  END DO
  ! ---------------------------------------------------------
  ! ---------------------------------------------------------
  ! ***   On repasse les S dans la variable qpr
  ! ***   On repasse les S dans la variable q directement 14/10/94

  DO l = 1, llm
    DO j = 1, jjp1
      DO i = 1, iip1
        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)
        q(i, j, llm+1-l, 4) = sxx(i, j, l)/sm(i, j, l)
        q(i, j, llm+1-l, 5) = sxy(i, j, l)/sm(i, j, l)
        q(i, j, llm+1-l, 6) = sxz(i, j, l)/sm(i, j, l)
        q(i, j, llm+1-l, 7) = syy(i, j, l)/sm(i, j, l)
        q(i, j, llm+1-l, 8) = syz(i, j, l)/sm(i, j, l)
        q(i, j, llm+1-l, 9) = szz(i, j, l)/sm(i, j, l)
      END DO
    END DO
  END DO

  ! ---------------------------------------------------------
  ! go to  777
  ! filtrages aux poles

  ! Traitements specifiques au pole

  ! filtrages aux poles
  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
    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 DO
  DO l = 1, llm
    DO i = 1, iim
      q(i, 1, llm+1-l, 4) = 0.
      q(i, jjp1, llm+1-l, 4) = 0.
      q(i, 1, llm+1-l, 5) = 0.
      q(i, jjp1, llm+1-l, 5) = 0.
      q(i, 1, llm+1-l, 6) = 0.
      q(i, jjp1, llm+1-l, 6) = 0.
      q(i, 1, llm+1-l, 7) = 0.
      q(i, jjp1, llm+1-l, 7) = 0.
      q(i, 1, llm+1-l, 8) = 0.
      q(i, jjp1, llm+1-l, 8) = 0.
      q(i, 1, llm+1-l, 9) = 0.
      q(i, jjp1, llm+1-l, 9) = 0.
    END DO
  END DO

  ! bouclage en longitude
  DO l = 1, llm
    DO j = 1, jjp1
      q(iip1, j, l, 0) = q(1, j, l, 0)
      q(iip1, j, llm+1-l, 0) = q(1, j, llm+1-l, 0)
      q(iip1, j, llm+1-l, 1) = q(1, j, llm+1-l, 1)
      q(iip1, j, llm+1-l, 2) = q(1, j, llm+1-l, 2)
      q(iip1, j, llm+1-l, 3) = q(1, j, llm+1-l, 3)
      q(iip1, j, llm+1-l, 4) = q(1, j, llm+1-l, 4)
      q(iip1, j, llm+1-l, 5) = q(1, j, llm+1-l, 5)
      q(iip1, j, llm+1-l, 6) = q(1, j, llm+1-l, 6)
      q(iip1, j, llm+1-l, 7) = q(1, j, llm+1-l, 7)
      q(iip1, j, llm+1-l, 8) = q(1, j, llm+1-l, 8)
      q(iip1, j, llm+1-l, 9) = q(1, j, llm+1-l, 9)
    END DO
  END DO
  DO l = 1, llm
    DO j = 2, jjm
      DO i = 1, iip1
        IF (q(i,j,l,0)<0.) THEN
          PRINT *, '------------ BIP-----------'
          PRINT *, 'S0(', i, j, l, ')=', q(i, j, l, 0), q(i, j-1, l, 0)
          PRINT *, 'SX(', i, j, l, ')=', q(i, j, l, 1)
          PRINT *, 'SY(', i, j, l, ')=', q(i, j, l, 2), q(i, j-1, l, 2)
          PRINT *, 'SZ(', i, j, l, ')=', q(i, j, l, 3)
          q(i, j, l, 0) = 0.
        END IF
      END DO
    END DO
    DO j = 1, jjp1, jjm
      DO i = 1, iip1
        IF (q(i,j,l,0)<0.) THEN
          PRINT *, '------------ BIP 2-----------'
          PRINT *, 'S0(', i, j, l, ')=', q(i, j, l, 0)
          PRINT *, 'SX(', i, j, l, ')=', q(i, j, l, 1)
          PRINT *, 'SY(', i, j, l, ')=', q(i, j, l, 2)
          PRINT *, 'SZ(', i, j, l, ')=', q(i, j, l, 3)

          q(i, j, l, 0) = 0.
          ! STOP
        END IF
      END DO
    END DO
  END DO
  RETURN
END SUBROUTINE prather
1	guez	3
2	guez	81	! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/prather.F,v 1.1.1.1 2004/05/19
3			! 12:53:07 lmdzadmin Exp $
4	guez	3
5	guez	81	SUBROUTINE prather(q, w, masse, pbaru, pbarv, nt, dt)
6	guez	139
7			USE comconst
8	guez	81	USE dimens_m
9			USE disvert_m
10	guez	139	USE dynetat0_m, only: rlonv, rlonu
11	guez	81	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 de prather
22			! Ref :
23	guez	3
24	guez	81	! ************************************************
25			! q,w,pext,pbaru et pbarv : arguments d'entree pour le s-pg
26	guez	3
27	guez	81	! =======================================================================
28	guez	3
29
30
31	guez	81	! Arguments:
32			! ----------
33			INTEGER iq, nt
34			REAL, INTENT (IN) :: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
35			REAL masse(iip1, jjp1, llm)
36			REAL q(iip1, jjp1, llm, 0:9)
37			REAL w(ip1jmp1, llm)
38			INTEGER ordre, ilim
39	guez	3
40	guez	81	! Local:
41			! ------
42			LOGICAL limit
43			REAL zq(iip1, jjp1, llm)
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 sxx(iip1, jjp1, llm)
48			REAL sxy(iip1, jjp1, llm)
49			REAL sxz(iip1, jjp1, llm)
50			REAL syy(iip1, jjp1, llm)
51			REAL syz(iip1, jjp1, llm)
52			REAL szz(iip1, jjp1, llm), zz
53			INTEGER i, j, l, indice
54			REAL sxn(iip1), sxs(iip1)
55	guez	3
56	guez	81	REAL sinlon(iip1), sinlondlon(iip1)
57			REAL coslon(iip1), coslondlon(iip1)
58			REAL qmin, qmax
59			SAVE qmin, qmax
60			SAVE sinlon, coslon, sinlondlon, coslondlon
61			REAL dyn1, dyn2, dys1, dys2, qpn, qps, dqzpn, dqzps
62			REAL masn, mass
63	guez	3
64	guez	81	REAL ssum
65			INTEGER ismax, ismin
66			EXTERNAL ssum, convflu, ismin, ismax
67			LOGICAL first
68			SAVE first
69			EXTERNAL advxp, advyp, advzp
70	guez	3
71
72	guez	81	DATA first/.TRUE./
73			DATA qmin, qmax/ -1.E33, 1.E33/
74	guez	3
75
76	guez	81	! ==========================================================================
77			! ==========================================================================
78			! MODIFICATION POUR PAS DE TEMPS ADAPTATIF, dtvr remplace par dt
79			! ==========================================================================
80			! ==========================================================================
81			REAL dt
82			! ==========================================================================
83			limit = .TRUE.
84	guez	3
85	guez	81	IF (first) THEN
86			PRINT *, 'SCHEMA PRATHER'
87			first = .FALSE.
88			DO i = 2, iip1
89			coslon(i) = cos(rlonv(i))
90			sinlon(i) = sin(rlonv(i))
91			coslondlon(i) = coslon(i)*(rlonu(i)-rlonu(i-1))/pi
92			sinlondlon(i) = sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
93			END DO
94			coslon(1) = coslon(iip1)
95			coslondlon(1) = coslondlon(iip1)
96			sinlon(1) = sinlon(iip1)
97			sinlondlon(1) = sinlondlon(iip1)
98	guez	3
99	guez	81	DO l = 1, llm
100			DO j = 1, jjp1
101			DO i = 1, iip1
102			q(i, j, l, 1) = 0.
103			q(i, j, l, 2) = 0.
104			q(i, j, l, 3) = 0.
105			q(i, j, l, 4) = 0.
106			q(i, j, l, 5) = 0.
107			q(i, j, l, 6) = 0.
108			q(i, j, l, 7) = 0.
109			q(i, j, l, 8) = 0.
110			q(i, j, l, 9) = 0.
111			END DO
112			END DO
113			END DO
114			END IF
115			! Fin modif Fred
116	guez	3
117	guez	81	! *** On calcule la masse d'air en kg
118	guez	3
119	guez	81	DO l = 1, llm
120			DO j = 1, jjp1
121			DO i = 1, iip1
122			sm(i, j, llm+1-l) = masse(i, j, l)
123			END DO
124			END DO
125			END DO
126	guez	3
127	guez	81	! *** q contient les qqtes de traceur avant l'advection
128	guez	3
129	guez	81	! *** Affectation des tableaux S a partir de Q
130
131			DO l = 1, llm
132			DO j = 1, jjp1
133			DO i = 1, iip1
134			s0(i, j, l) = q(i, j, llm+1-l, 0)*sm(i, j, l)
135			sx(i, j, l) = q(i, j, llm+1-l, 1)*sm(i, j, l)
136			sy(i, j, l) = q(i, j, llm+1-l, 2)*sm(i, j, l)
137			sz(i, j, l) = q(i, j, llm+1-l, 3)*sm(i, j, l)
138			sxx(i, j, l) = q(i, j, llm+1-l, 4)*sm(i, j, l)
139			sxy(i, j, l) = q(i, j, llm+1-l, 5)*sm(i, j, l)
140			sxz(i, j, l) = q(i, j, llm+1-l, 6)*sm(i, j, l)
141			syy(i, j, l) = q(i, j, llm+1-l, 7)*sm(i, j, l)
142			syz(i, j, l) = q(i, j, llm+1-l, 8)*sm(i, j, l)
143			szz(i, j, l) = q(i, j, llm+1-l, 9)*sm(i, j, l)
144	guez	3	END DO
145	guez	81	END DO
146			END DO
147			! *** Appel des subroutines d'advection en X, en Y et en Z
148			! *** Advection avec "time-splitting"
149	guez	3
150	guez	81	! -----------------------------------------------------------
151			DO indice = 1, nt
152			CALL advxp(limit, 0.5*dt, pbaru, sm, s0, sx, sy, sz, sxx, sxy, sxz, syy, &
153			syz, szz, 1)
154			END DO
155			DO l = 1, llm
156			DO i = 1, iip1
157			sy(i, 1, l) = 0.
158			sy(i, jjp1, l) = 0.
159			END DO
160			END DO
161			! ---------------------------------------------------------
162			CALL advyp(limit, .5dtnt, pbarv, sm, s0, sx, sy, sz, sxx, sxy, sxz, syy, &
163			syz, szz, 1)
164			! ---------------------------------------------------------
165	guez	3
166	guez	81	! ---------------------------------------------------------
167			DO j = 1, jjp1
168			DO i = 1, iip1
169			sz(i, j, 1) = 0.
170			sz(i, j, llm) = 0.
171			sxz(i, j, 1) = 0.
172			sxz(i, j, llm) = 0.
173			syz(i, j, 1) = 0.
174			syz(i, j, llm) = 0.
175			szz(i, j, 1) = 0.
176			szz(i, j, llm) = 0.
177			END DO
178			END DO
179			CALL advzp(limit, dt*nt, w, sm, s0, sx, sy, sz, sxx, sxy, sxz, syy, syz, &
180			szz, 1)
181			DO l = 1, llm
182			DO i = 1, iip1
183			sy(i, 1, l) = 0.
184			sy(i, jjp1, l) = 0.
185			END DO
186			END DO
187
188			! ---------------------------------------------------------
189
190			! ---------------------------------------------------------
191			CALL advyp(limit, .5dtnt, pbarv, sm, s0, sx, sy, sz, sxx, sxy, sxz, syy, &
192			syz, szz, 1)
193			! ---------------------------------------------------------
194			DO l = 1, llm
195			DO j = 1, jjp1
196			s0(iip1, j, l) = s0(1, j, l)
197			sx(iip1, j, l) = sx(1, j, l)
198			sy(iip1, j, l) = sy(1, j, l)
199			sz(iip1, j, l) = sz(1, j, l)
200			sxx(iip1, j, l) = sxx(1, j, l)
201			sxy(iip1, j, l) = sxy(1, j, l)
202			sxz(iip1, j, l) = sxz(1, j, l)
203			syy(iip1, j, l) = syy(1, j, l)
204			syz(iip1, j, l) = syz(1, j, l)
205			szz(iip1, j, l) = szz(1, j, l)
206			END DO
207			END DO
208			DO indice = 1, nt
209			CALL advxp(limit, 0.5*dt, pbaru, sm, s0, sx, sy, sz, sxx, sxy, sxz, syy, &
210			syz, szz, 1)
211			END DO
212			! ---------------------------------------------------------
213			! ---------------------------------------------------------
214			! *** On repasse les S dans la variable qpr
215			! *** On repasse les S dans la variable q directement 14/10/94
216
217			DO l = 1, llm
218			DO j = 1, jjp1
219			DO i = 1, iip1
220			q(i, j, llm+1-l, 0) = s0(i, j, l)/sm(i, j, l)
221			q(i, j, llm+1-l, 1) = sx(i, j, l)/sm(i, j, l)
222			q(i, j, llm+1-l, 2) = sy(i, j, l)/sm(i, j, l)
223			q(i, j, llm+1-l, 3) = sz(i, j, l)/sm(i, j, l)
224			q(i, j, llm+1-l, 4) = sxx(i, j, l)/sm(i, j, l)
225			q(i, j, llm+1-l, 5) = sxy(i, j, l)/sm(i, j, l)
226			q(i, j, llm+1-l, 6) = sxz(i, j, l)/sm(i, j, l)
227			q(i, j, llm+1-l, 7) = syy(i, j, l)/sm(i, j, l)
228			q(i, j, llm+1-l, 8) = syz(i, j, l)/sm(i, j, l)
229			q(i, j, llm+1-l, 9) = szz(i, j, l)/sm(i, j, l)
230			END DO
231			END DO
232			END DO
233
234			! ---------------------------------------------------------
235			! go to 777
236			! filtrages aux poles
237
238			! Traitements specifiques au pole
239
240			! filtrages aux poles
241			DO l = 1, llm
242			! filtrages aux poles
243			masn = ssum(iim, sm(1,1,l), 1)
244			mass = ssum(iim, sm(1,jjp1,l), 1)
245			qpn = ssum(iim, s0(1,1,l), 1)/masn
246			qps = ssum(iim, s0(1,jjp1,l), 1)/mass
247			dqzpn = ssum(iim, sz(1,1,l), 1)/masn
248			dqzps = ssum(iim, sz(1,jjp1,l), 1)/mass
249			DO i = 1, iip1
250			q(i, 1, llm+1-l, 3) = dqzpn
251			q(i, jjp1, llm+1-l, 3) = dqzps
252			q(i, 1, llm+1-l, 0) = qpn
253			q(i, jjp1, llm+1-l, 0) = qps
254			END DO
255			dyn1 = 0.
256			dys1 = 0.
257			dyn2 = 0.
258			dys2 = 0.
259			DO i = 1, iim
260			zz = s0(i, 2, l)/sm(i, 2, l) - q(i, 1, llm+1-l, 0)
261			dyn1 = dyn1 + sinlondlon(i)*zz
262			dyn2 = dyn2 + coslondlon(i)*zz
263			zz = q(i, jjp1, llm+1-l, 0) - s0(i, jjm, l)/sm(i, jjm, l)
264			dys1 = dys1 + sinlondlon(i)*zz
265			dys2 = dys2 + coslondlon(i)*zz
266			END DO
267			DO i = 1, iim
268			q(i, 1, llm+1-l, 2) = (sinlon(i)dyn1+coslon(i)dyn2)/2.
269			q(i, 1, llm+1-l, 0) = q(i, 1, llm+1-l, 0) + q(i, 1, llm+1-l, 2)
270			q(i, jjp1, llm+1-l, 2) = (sinlon(i)dys1+coslon(i)dys2)/2.
271			q(i, jjp1, llm+1-l, 0) = q(i, jjp1, llm+1-l, 0) - &
272			q(i, jjp1, llm+1-l, 2)
273			END DO
274			q(iip1, 1, llm+1-l, 0) = q(1, 1, llm+1-l, 0)
275			q(iip1, jjp1, llm+1-l, 0) = q(1, jjp1, llm+1-l, 0)
276			DO i = 1, iim
277			sxn(i) = q(i+1, 1, llm+1-l, 0) - q(i, 1, llm+1-l, 0)
278			sxs(i) = q(i+1, jjp1, llm+1-l, 0) - q(i, jjp1, llm+1-l, 0)
279			END DO
280			sxn(iip1) = sxn(1)
281			sxs(iip1) = sxs(1)
282			DO i = 1, iim
283			q(i+1, 1, llm+1-l, 1) = 0.25*(sxn(i)+sxn(i+1))
284			q(i+1, jjp1, llm+1-l, 1) = 0.25*(sxs(i)+sxs(i+1))
285			END DO
286			q(1, 1, llm+1-l, 1) = q(iip1, 1, llm+1-l, 1)
287			q(1, jjp1, llm+1-l, 1) = q(iip1, jjp1, llm+1-l, 1)
288			END DO
289			DO l = 1, llm
290			DO i = 1, iim
291			q(i, 1, llm+1-l, 4) = 0.
292			q(i, jjp1, llm+1-l, 4) = 0.
293			q(i, 1, llm+1-l, 5) = 0.
294			q(i, jjp1, llm+1-l, 5) = 0.
295			q(i, 1, llm+1-l, 6) = 0.
296			q(i, jjp1, llm+1-l, 6) = 0.
297			q(i, 1, llm+1-l, 7) = 0.
298			q(i, jjp1, llm+1-l, 7) = 0.
299			q(i, 1, llm+1-l, 8) = 0.
300			q(i, jjp1, llm+1-l, 8) = 0.
301			q(i, 1, llm+1-l, 9) = 0.
302			q(i, jjp1, llm+1-l, 9) = 0.
303			END DO
304			END DO
305
306			! bouclage en longitude
307			DO l = 1, llm
308			DO j = 1, jjp1
309			q(iip1, j, l, 0) = q(1, j, l, 0)
310			q(iip1, j, llm+1-l, 0) = q(1, j, llm+1-l, 0)
311			q(iip1, j, llm+1-l, 1) = q(1, j, llm+1-l, 1)
312			q(iip1, j, llm+1-l, 2) = q(1, j, llm+1-l, 2)
313			q(iip1, j, llm+1-l, 3) = q(1, j, llm+1-l, 3)
314			q(iip1, j, llm+1-l, 4) = q(1, j, llm+1-l, 4)
315			q(iip1, j, llm+1-l, 5) = q(1, j, llm+1-l, 5)
316			q(iip1, j, llm+1-l, 6) = q(1, j, llm+1-l, 6)
317			q(iip1, j, llm+1-l, 7) = q(1, j, llm+1-l, 7)
318			q(iip1, j, llm+1-l, 8) = q(1, j, llm+1-l, 8)
319			q(iip1, j, llm+1-l, 9) = q(1, j, llm+1-l, 9)
320			END DO
321			END DO
322			DO l = 1, llm
323			DO j = 2, jjm
324			DO i = 1, iip1
325			IF (q(i,j,l,0)<0.) THEN
326			PRINT *, '------------ BIP-----------'
327			PRINT *, 'S0(', i, j, l, ')=', q(i, j, l, 0), q(i, j-1, l, 0)
328			PRINT *, 'SX(', i, j, l, ')=', q(i, j, l, 1)
329			PRINT *, 'SY(', i, j, l, ')=', q(i, j, l, 2), q(i, j-1, l, 2)
330			PRINT *, 'SZ(', i, j, l, ')=', q(i, j, l, 3)
331			q(i, j, l, 0) = 0.
332			END IF
333			END DO
334			END DO
335			DO j = 1, jjp1, jjm
336			DO i = 1, iip1
337			IF (q(i,j,l,0)<0.) THEN
338			PRINT *, '------------ BIP 2-----------'
339			PRINT *, 'S0(', i, j, l, ')=', q(i, j, l, 0)
340			PRINT *, 'SX(', i, j, l, ')=', q(i, j, l, 1)
341			PRINT *, 'SY(', i, j, l, ')=', q(i, j, l, 2)
342			PRINT *, 'SZ(', i, j, l, ')=', q(i, j, l, 3)
343
344			q(i, j, l, 0) = 0.
345			! STOP
346			END IF
347			END DO
348			END DO
349			END DO
350			RETURN
351			END SUBROUTINE prather