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
2	! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/prather.F,v 1.1.1.1 2004/05/19
3	! 12:53:07 lmdzadmin Exp $
4
5	SUBROUTINE prather(q, w, masse, pbaru, pbarv, nt, dt)
6
7	USE comconst
8	USE dimens_m
9	USE disvert_m
10	USE dynetat0_m, only: rlonv, rlonu
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 de prather
22	! Ref :
23
24	! ************************************************
25	! q,w,pext,pbaru et pbarv : arguments d'entree pour le s-pg
26
27	! =======================================================================
28
29
30
31	! 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
40	! 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
56	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
64	REAL ssum
65	INTEGER ismax, ismin
66	EXTERNAL ssum, convflu, ismin, ismax
67	LOGICAL first
68	SAVE first
69	EXTERNAL advxp, advyp, advzp
70
71
72	DATA first/.TRUE./
73	DATA qmin, qmax/ -1.E33, 1.E33/
74
75
76	! ==========================================================================
77	! ==========================================================================
78	! MODIFICATION POUR PAS DE TEMPS ADAPTATIF, dtvr remplace par dt
79	! ==========================================================================
80	! ==========================================================================
81	REAL dt
82	! ==========================================================================
83	limit = .TRUE.
84
85	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
99	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
117	! *** On calcule la masse d'air en kg
118
119	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
127	! *** q contient les qqtes de traceur avant l'advection
128
129	! *** 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	END DO
145	END DO
146	END DO
147	! *** Appel des subroutines d'advection en X, en Y et en Z
148	! *** Advection avec "time-splitting"
149
150	! -----------------------------------------------------------
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
166	! ---------------------------------------------------------
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