trunk/dyn3d/prather.f90


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