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 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)

  ! 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 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)
  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./

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