libf/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 dimens_m
      use paramet_m
      use comconst
      use comvert
      use comgeom
      IMPLICIT NONE

c=======================================================================
c   Adaptation LMDZ:  A.Armengaud (LGGE)
c   ----------------
c
c   ************************************************
c   Transport des traceurs par la methode de prather
c   Ref : 
c
c   ************************************************
c   q,w,pext,pbaru et pbarv : arguments d'entree  pour le s-pg
c
c=======================================================================


c   Arguments:
c   ----------
      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

c   Local:
c   ------
      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
c
      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/


c==========================================================================
c==========================================================================
c     MODIFICATION POUR PAS DE TEMPS ADAPTATIF, dtvr remplace par dt
c==========================================================================
c==========================================================================
      REAL dt
c==========================================================================
      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
         enddo
         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.
        ENDDO
        ENDDO
        ENDDO
      endif
c   Fin modif Fred

c *** 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)
         ENDDO
        ENDDO
       ENDDO

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

c *** 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)
         ENDDO
        ENDDO
       ENDDO
c *** Appel des subroutines d'advection en X, en Y et en Z
c *** Advection avec "time-splitting"
      
c-----------------------------------------------------------
       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.
        enddo
        enddo
c---------------------------------------------------------
       call advyp( limit,.5*dt*nt,pbarv,sm,s0,sx,sy,sz
     .             ,sxx,sxy,sxz,syy,syz,szz,1 )
c---------------------------------------------------------

c---------------------------------------------------------
       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.
          enddo
       enddo
       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.
        enddo
        enddo

c---------------------------------------------------------

c---------------------------------------------------------
       call advyp( limit,.5*dt*nt,pbarv,sm,s0,sx,sy,sz
     .             ,sxx,sxy,sxz,syy,syz,szz,1 )
c---------------------------------------------------------
       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 )
        ENDDO
       ENDDO
       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
c---------------------------------------------------------
c---------------------------------------------------------
c ***   On repasse les S dans la variable qpr
c ***   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)
      ENDDO
      ENDDO
      ENDDO

c---------------------------------------------------------
c      go to  777
c   filtrages aux poles

c Traitements specifiques au pole

c   filtrages aux poles
         DO l=1,llm
c   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
         enddo
c       enddo
c         print*,'qpn',qpn,'qps',qps
c          print*,'dqzpn',dqzpn,'dqzps',dqzps
c       enddo
           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
        enddo
         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)
         enddo
      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)
      enddo
      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)
        enddo
         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.
          enddo
         ENDDO

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