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
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	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
14	! =======================================================================
15	! Adaptation LMDZ: A.Armengaud (LGGE)
16	! ----------------
17
18	! ************************************************
19	! Transport des traceurs par la methode de prather
20	! Ref :
21
22	! ************************************************
23	! q,w,pext,pbaru et pbarv : arguments d'entree pour le s-pg
24
25	! =======================================================================
26
27
28
29	! 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
38	! 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
54	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
62	REAL ssum
63	INTEGER ismax, ismin
64	EXTERNAL ssum, convflu, ismin, ismax
65	LOGICAL first
66	SAVE first
67	EXTERNAL advxp, advyp, advzp
68
69
70	DATA first/.TRUE./
71	DATA qmin, qmax/ -1.E33, 1.E33/
72
73
74	! ==========================================================================
75	! ==========================================================================
76	! MODIFICATION POUR PAS DE TEMPS ADAPTATIF, dtvr remplace par dt
77	! ==========================================================================
78	! ==========================================================================
79	REAL dt
80	! ==========================================================================
81	limit = .TRUE.
82
83	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
97	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
115	! *** On calcule la masse d'air en kg
116
117	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
125	! *** q contient les qqtes de traceur avant l'advection
126
127	! *** 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	END DO
143	END DO
144	END DO
145	! *** Appel des subroutines d'advection en X, en Y et en Z
146	! *** Advection avec "time-splitting"
147
148	! -----------------------------------------------------------
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
164	! ---------------------------------------------------------
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