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