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	!
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	REAL, intent(in):: pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm )
31	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