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
           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)
                     q(i,j,l,0)=0.
               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			dyn1=0.
252			dys1=0.
253			dyn2=0.
254			dys2=0.
255			do i=1,iim
256			zz=s0(i,2,l)/sm(i,2,l)-q(i,1,llm+1-l,0)
257			dyn1=dyn1+sinlondlon(i)*zz
258			dyn2=dyn2+coslondlon(i)*zz
259			zz=q(i,jjp1,llm+1-l,0)-s0(i,jjm,l)/sm(i,jjm,l)
260			dys1=dys1+sinlondlon(i)*zz
261			dys2=dys2+coslondlon(i)*zz
262			enddo
263			do i=1,iim
264			q(i,1,llm+1-l,2)=
265			$ (sinlon(i)dyn1+coslon(i)dyn2)/2.
266			q(i,1,llm+1-l,0)=q(i,1,llm+1-l,0)
267			$ +q(i,1,llm+1-l,2)
268			q(i,jjp1,llm+1-l,2)=
269			$ (sinlon(i)dys1+coslon(i)dys2)/2.
270			q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0)
271			$ -q(i,jjp1,llm+1-l,2)
272			enddo
273			q(iip1,1,llm+1-l,0)=q(1,1,llm+1-l,0)
274			q(iip1,jjp1,llm+1-l,0)=q(1,jjp1,llm+1-l,0)
275			do i=1,iim
276			sxn(i)=q(i+1,1,llm+1-l,0)-q(i,1,llm+1-l,0)
277			sxs(i)=q(i+1,jjp1,llm+1-l,0)-q(i,jjp1,llm+1-l,0)
278			enddo
279			sxn(iip1)=sxn(1)
280			sxs(iip1)=sxs(1)
281			do i=1,iim
282			q(i+1,1,llm+1-l,1)=0.25*(sxn(i)+sxn(i+1))
283			q(i+1,jjp1,llm+1-l,1)=0.25*(sxs(i)+sxs(i+1))
284			END DO
285			q(1,1,llm+1-l,1)=q(iip1,1,llm+1-l,1)
286			q(1,jjp1,llm+1-l,1)=
287			$ q(iip1,jjp1,llm+1-l,1)
288			enddo
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			enddo
304			ENDDO
305
306			777 continue
307			c
308			c bouclage en longitude
309			do l=1,llm
310			do j=1,jjp1
311			q(iip1,j,l,0)=q(1,j,l,0)
312			q(iip1,j,llm+1-l,0)=q(1,j,llm+1-l,0)
313			q(iip1,j,llm+1-l,1)=q(1,j,llm+1-l,1)
314			q(iip1,j,llm+1-l,2)=q(1,j,llm+1-l,2)
315			q(iip1,j,llm+1-l,3)=q(1,j,llm+1-l,3)
316			q(iip1,j,llm+1-l,4)=q(1,j,llm+1-l,4)
317			q(iip1,j,llm+1-l,5)=q(1,j,llm+1-l,5)
318			q(iip1,j,llm+1-l,6)=q(1,j,llm+1-l,6)
319			q(iip1,j,llm+1-l,7)=q(1,j,llm+1-l,7)
320			q(iip1,j,llm+1-l,8)=q(1,j,llm+1-l,8)
321			q(iip1,j,llm+1-l,9)=q(1,j,llm+1-l,9)
322			enddo
323			enddo
324			DO l = 1,llm
325			DO j = 2,jjm
326			DO i = 1,iip1
327			IF (q(i,j,l,0).lt.0.) THEN
328			PRINT*,'------------ BIP-----------'
329			PRINT*,'S0(',i,j,l,')=',q(i,j,l,0),
330			$ q(i,j-1,l,0)
331			PRINT*,'SX(',i,j,l,')=',q(i,j,l,1)
332			PRINT*,'SY(',i,j,l,')=',q(i,j,l,2),
333			$ q(i,j-1,l,2)
334			PRINT*,'SZ(',i,j,l,')=',q(i,j,l,3)
335			q(i,j,l,0)=0.
336			ENDIF
337			ENDDO
338			ENDDO
339			do j=1,jjp1,jjm
340			do i=1,iip1
341			IF (q(i,j,l,0).lt.0.) THEN
342			PRINT*,'------------ BIP 2-----------'
343			PRINT*,'S0(',i,j,l,')=',q(i,j,l,0)
344			PRINT*,'SX(',i,j,l,')=',q(i,j,l,1)
345			PRINT*,'SY(',i,j,l,')=',q(i,j,l,2)
346			PRINT*,'SZ(',i,j,l,')=',q(i,j,l,3)
347
348			q(i,j,l,0)=0.
349			c STOP
350			ENDIF
351			enddo
352			enddo
353			ENDDO
354			RETURN
355			END