dyn3d/Vlsplt/vly.f

      SUBROUTINE vly(q,pente_max,masse,masse_adv_v)
!
!     Auteurs:   P.Le Van, F.Hourdin, F.Forget
!
!    ***********************************************************
!     Shema  d'advection " pseudo amont " .
!    *************************************************************
!     q,masse_adv_v,w sont des arguments d'entree  pour le s-pg ....
!     dq              sont des arguments de sortie pour le s-pg ....
!
!
!   ----------------------------------------------------------------
      use dimens_m
      use paramet_m
      use comconst
      use comvert
      use logic
      use comgeom
      IMPLICIT NONE
!
!
!
!   Arguments:
!   ----------
      REAL masse(ip1jmp1,llm),pente_max
      REAL masse_adv_v( ip1jm,llm)
      REAL q(ip1jmp1,llm), dq( ip1jmp1,llm)
!
!      Local
!   ---------
!
      INTEGER i,ij,l
!
      REAL airej2,airejjm,airescb(iim),airesch(iim)
      REAL dyq(ip1jmp1,llm),dyqv(ip1jm),zdvm(ip1jmp1,llm)
      REAL adyqv(ip1jm),dyqmax(ip1jmp1)
      REAL qbyv(ip1jm,llm)

      REAL qpns,qpsn,apn,aps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
!     REAL newq,oldmasse
      Logical extremum,first,testcpu
      REAL temps0,temps1,temps2,temps3,temps4,temps5,second
      SAVE temps0,temps1,temps2,temps3,temps4,temps5
      SAVE first,testcpu

      REAL convpn,convps,convmpn,convmps
      real massepn,masseps,qpn,qps
      REAL sinlon(iip1),sinlondlon(iip1)
      REAL coslon(iip1),coslondlon(iip1)
      SAVE sinlon,coslon,sinlondlon,coslondlon
      SAVE airej2,airejjm
!
!
      REAL      SSUM

      DATA first,testcpu/.true.,.false./
      DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./

      IF(first) THEN
         PRINT*,'Shema  Amont nouveau  appele dans  Vanleer   '
         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)
         airej2 = SSUM( iim, aire(iip2), 1 )
         airejjm= SSUM( iim, aire(ip1jm -iim), 1 )
      ENDIF

!
!PRINT*,'CALCUL EN LATITUDE'

      DO l = 1, llm
!
!   --------------------------------
!      CALCUL EN LATITUDE
!   --------------------------------

!   On commence par calculer la valeur du traceur moyenne sur le
!     premier cercle
!   de latitude autour du pole (qpns pour le pole nord et qpsn pour
!    le pole nord) qui sera utilisee pour evaluer les pentes au pole.

      DO i = 1, iim
      airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
      airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
      ENDDO
      qpns   = SSUM( iim,  airescb ,1 ) / airej2
      qpsn   = SSUM( iim,  airesch ,1 ) / airejjm

!   calcul des pentes aux points v

      DO ij=1,ip1jm
         dyqv(ij)=q(ij,l)-q(ij+iip1,l)
         adyqv(ij)=abs(dyqv(ij))
      ENDDO

!   calcul des pentes aux points scalaires

      DO ij=iip2,ip1jm
         dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
         dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
         dyqmax(ij)=pente_max*dyqmax(ij)
      ENDDO

!   calcul des pentes aux poles

      DO ij=1,iip1
         dyq(ij,l)=qpns-q(ij+iip1,l)
         dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
      ENDDO

!   filtrage de la derivee
      dyn1=0.
      dys1=0.
      dyn2=0.
      dys2=0.
      DO ij=1,iim
         dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
         dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
         dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
         dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
      ENDDO
      DO ij=1,iip1
         dyq(ij,l)=dyn1*sinlon(ij)+dyn2*coslon(ij)
         dyq(ip1jm+ij,l)=dys1*sinlon(ij)+dys2*coslon(ij)
      ENDDO

!   calcul des pentes limites aux poles

      goto 8888
      fn=1.
      fs=1.
      DO ij=1,iim
         IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
            fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
         ENDIF
      IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
         fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
         ENDIF
      ENDDO
      DO ij=1,iip1
         dyq(ij,l)=fn*dyq(ij,l)
         dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
      ENDDO
 8888   continue
      DO ij=1,iip1
         dyq(ij,l)=0.
         dyq(ip1jm+ij,l)=0.
      ENDDO

!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!  En memoire de dIFferents tests sur la
!  limitation des pentes aux poles.
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!     PRINT*,dyq(1)
!     PRINT*,dyqv(iip1+1)
!     apn=abs(dyq(1)/dyqv(iip1+1))
!     PRINT*,dyq(ip1jm+1)
!     PRINT*,dyqv(ip1jm-iip1+1)
!     aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
!     DO ij=2,iim
!        apn=amax1(abs(dyq(ij)/dyqv(ij)),apn)
!        aps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps)
!     ENDDO
!     apn=min(pente_max/apn,1.)
!     aps=min(pente_max/aps,1.)
!
!
!   cas ou on a un extremum au pole
!
!     IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
!    &   apn=0.
!     IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
!    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
!    &   aps=0.
!
!   limitation des pentes aux poles
!     DO ij=1,iip1
!        dyq(ij)=apn*dyq(ij)
!        dyq(ip1jm+ij)=aps*dyq(ip1jm+ij)
!     ENDDO
!
!   test
!      DO ij=1,iip1
!         dyq(iip1+ij)=0.
!         dyq(ip1jm+ij-iip1)=0.
!      ENDDO
!      DO ij=1,ip1jmp1
!         dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
!      ENDDO
!
! changement 10 07 96
!     IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
!    &   THEN
!        DO ij=1,iip1
!           dyqmax(ij)=0.
!        ENDDO
!     ELSE
!        DO ij=1,iip1
!           dyqmax(ij)=pente_max*abs(dyqv(ij))
!        ENDDO
!     ENDIF
!
!     IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
!    & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
!    &THEN
!        DO ij=ip1jm+1,ip1jmp1
!           dyqmax(ij)=0.
!        ENDDO
!     ELSE
!        DO ij=ip1jm+1,ip1jmp1
!           dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
!        ENDDO
!     ENDIF
!   fin changement 10 07 96
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

!   calcul des pentes limitees

      DO ij=iip2,ip1jm
         IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
            dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
         ELSE
            dyq(ij,l)=0.
         ENDIF
      ENDDO

      ENDDO

      DO l=1,llm
       DO ij=1,ip1jm
          IF(masse_adv_v(ij,l).gt.0) THEN
              qbyv(ij,l)=q(ij+iip1,l)+dyq(ij+iip1,l)*                   &
     &                   0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l))
          ELSE
              qbyv(ij,l)=q(ij,l)-dyq(ij,l)*                             &
     &                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l))
          ENDIF
          qbyv(ij,l)=masse_adv_v(ij,l)*qbyv(ij,l)
       ENDDO
      ENDDO


      DO l=1,llm
         DO ij=iip2,ip1jm
            newmasse=masse(ij,l)                                        &
     &      +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
            q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))    &
     &         /newmasse
            masse(ij,l)=newmasse
         ENDDO
!.-. ancienne version
!        convpn=SSUM(iim,qbyv(1,l),1)/apoln
!        convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln

         convpn=SSUM(iim,qbyv(1,l),1)
         convmpn=ssum(iim,masse_adv_v(1,l),1)
         massepn=ssum(iim,masse(1,l),1)
         qpn=0.
         do ij=1,iim
            qpn=qpn+masse(ij,l)*q(ij,l)
         enddo
         qpn=(qpn+convpn)/(massepn+convmpn)
         do ij=1,iip1
            q(ij,l)=qpn
         enddo

!        convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
!        convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)/apols

         convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
         convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
         masseps=ssum(iim, masse(ip1jm+1,l),1)
         qps=0.
         do ij = ip1jm+1,ip1jmp1-1
            qps=qps+masse(ij,l)*q(ij,l)
         enddo
         qps=(qps+convps)/(masseps+convmps)
         do ij=ip1jm+1,ip1jmp1
            q(ij,l)=qps
         enddo

!.-. fin ancienne version

!._. nouvelle version
!        convpn=SSUM(iim,qbyv(1,l),1)
!        convmpn=ssum(iim,masse_adv_v(1,l),1)
!        oldmasse=ssum(iim,masse(1,l),1)
!        newmasse=oldmasse+convmpn
!        newq=(q(1,l)*oldmasse+convpn)/newmasse
!        newmasse=newmasse/apoln
!        DO ij = 1,iip1
!           q(ij,l)=newq
!           masse(ij,l)=newmasse*aire(ij)
!        ENDDO
!        convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
!        convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
!        oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
!        newmasse=oldmasse+convmps
!        newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
!        newmasse=newmasse/apols
!        DO ij = ip1jm+1,ip1jmp1
!           q(ij,l)=newq
!           masse(ij,l)=newmasse*aire(ij)
!        ENDDO
!._. fin nouvelle version
      ENDDO

      RETURN
      END
1	SUBROUTINE vly(q,pente_max,masse,masse_adv_v)
2	!
3	! Auteurs: P.Le Van, F.Hourdin, F.Forget
4	!
5	! ***********************************************************
6	! Shema d'advection " pseudo amont " .
7	! *************************************************************
8	! q,masse_adv_v,w sont des arguments d'entree pour le s-pg ....
9	! dq sont des arguments de sortie pour le s-pg ....
10	!
11	!
12	! ----------------------------------------------------------------
13	use dimens_m
14	use paramet_m
15	use comconst
16	use comvert
17	use logic
18	use comgeom
19	IMPLICIT NONE
20	!
21	!
22	!
23	! Arguments:
24	! ----------
25	REAL masse(ip1jmp1,llm),pente_max
26	REAL masse_adv_v( ip1jm,llm)
27	REAL q(ip1jmp1,llm), dq( ip1jmp1,llm)
28	!
29	! Local
30	! ---------
31	!
32	INTEGER i,ij,l
33	!
34	REAL airej2,airejjm,airescb(iim),airesch(iim)
35	REAL dyq(ip1jmp1,llm),dyqv(ip1jm),zdvm(ip1jmp1,llm)
36	REAL adyqv(ip1jm),dyqmax(ip1jmp1)
37	REAL qbyv(ip1jm,llm)
38
39	REAL qpns,qpsn,apn,aps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
40	! REAL newq,oldmasse
41	Logical extremum,first,testcpu
42	REAL temps0,temps1,temps2,temps3,temps4,temps5,second
43	SAVE temps0,temps1,temps2,temps3,temps4,temps5
44	SAVE first,testcpu
45
46	REAL convpn,convps,convmpn,convmps
47	real massepn,masseps,qpn,qps
48	REAL sinlon(iip1),sinlondlon(iip1)
49	REAL coslon(iip1),coslondlon(iip1)
50	SAVE sinlon,coslon,sinlondlon,coslondlon
51	SAVE airej2,airejjm
52	!
53	!
54	REAL SSUM
55
56	DATA first,testcpu/.true.,.false./
57	DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
58
59	IF(first) THEN
60	PRINT*,'Shema Amont nouveau appele dans Vanleer '
61	first=.false.
62	do i=2,iip1
63	coslon(i)=cos(rlonv(i))
64	sinlon(i)=sin(rlonv(i))
65	coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
66	sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
67	ENDDO
68	coslon(1)=coslon(iip1)
69	coslondlon(1)=coslondlon(iip1)
70	sinlon(1)=sinlon(iip1)
71	sinlondlon(1)=sinlondlon(iip1)
72	airej2 = SSUM( iim, aire(iip2), 1 )
73	airejjm= SSUM( iim, aire(ip1jm -iim), 1 )
74	ENDIF
75
76	!
77	!PRINT*,'CALCUL EN LATITUDE'
78
79	DO l = 1, llm
80	!
81	! --------------------------------
82	! CALCUL EN LATITUDE
83	! --------------------------------
84
85	! On commence par calculer la valeur du traceur moyenne sur le
86	! premier cercle
87	! de latitude autour du pole (qpns pour le pole nord et qpsn pour
88	! le pole nord) qui sera utilisee pour evaluer les pentes au pole.
89
90	DO i = 1, iim
91	airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
92	airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
93	ENDDO
94	qpns = SSUM( iim, airescb ,1 ) / airej2
95	qpsn = SSUM( iim, airesch ,1 ) / airejjm
96
97	! calcul des pentes aux points v
98
99	DO ij=1,ip1jm
100	dyqv(ij)=q(ij,l)-q(ij+iip1,l)
101	adyqv(ij)=abs(dyqv(ij))
102	ENDDO
103
104	! calcul des pentes aux points scalaires
105
106	DO ij=iip2,ip1jm
107	dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
108	dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
109	dyqmax(ij)=pente_max*dyqmax(ij)
110	ENDDO
111
112	! calcul des pentes aux poles
113
114	DO ij=1,iip1
115	dyq(ij,l)=qpns-q(ij+iip1,l)
116	dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
117	ENDDO
118
119	! filtrage de la derivee
120	dyn1=0.
121	dys1=0.
122	dyn2=0.
123	dys2=0.
124	DO ij=1,iim
125	dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
126	dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
127	dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
128	dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
129	ENDDO
130	DO ij=1,iip1
131	dyq(ij,l)=dyn1sinlon(ij)+dyn2coslon(ij)
132	dyq(ip1jm+ij,l)=dys1sinlon(ij)+dys2coslon(ij)
133	ENDDO
134
135	! calcul des pentes limites aux poles
136
137	goto 8888
138	fn=1.
139	fs=1.
140	DO ij=1,iim
141	IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
142	fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
143	ENDIF
144	IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
145	fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
146	ENDIF
147	ENDDO
148	DO ij=1,iip1
149	dyq(ij,l)=fn*dyq(ij,l)
150	dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
151	ENDDO
152	8888 continue
153	DO ij=1,iip1
154	dyq(ij,l)=0.
155	dyq(ip1jm+ij,l)=0.
156	ENDDO
157
158	!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
159	! En memoire de dIFferents tests sur la
160	! limitation des pentes aux poles.
161	!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
162	! PRINT*,dyq(1)
163	! PRINT*,dyqv(iip1+1)
164	! apn=abs(dyq(1)/dyqv(iip1+1))
165	! PRINT*,dyq(ip1jm+1)
166	! PRINT*,dyqv(ip1jm-iip1+1)
167	! aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
168	! DO ij=2,iim
169	! apn=amax1(abs(dyq(ij)/dyqv(ij)),apn)
170	! aps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps)
171	! ENDDO
172	! apn=min(pente_max/apn,1.)
173	! aps=min(pente_max/aps,1.)
174	!
175	!
176	! cas ou on a un extremum au pole
177	!
178	! IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
179	! & apn=0.
180	! IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
181	! & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
182	! & aps=0.
183	!
184	! limitation des pentes aux poles
185	! DO ij=1,iip1
186	! dyq(ij)=apn*dyq(ij)
187	! dyq(ip1jm+ij)=aps*dyq(ip1jm+ij)
188	! ENDDO
189	!
190	! test
191	! DO ij=1,iip1
192	! dyq(iip1+ij)=0.
193	! dyq(ip1jm+ij-iip1)=0.
194	! ENDDO
195	! DO ij=1,ip1jmp1
196	! dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
197	! ENDDO
198	!
199	! changement 10 07 96
200	! IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
201	! & THEN
202	! DO ij=1,iip1
203	! dyqmax(ij)=0.
204	! ENDDO
205	! ELSE
206	! DO ij=1,iip1
207	! dyqmax(ij)=pente_max*abs(dyqv(ij))
208	! ENDDO
209	! ENDIF
210	!
211	! IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
212	! & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
213	! &THEN
214	! DO ij=ip1jm+1,ip1jmp1
215	! dyqmax(ij)=0.
216	! ENDDO
217	! ELSE
218	! DO ij=ip1jm+1,ip1jmp1
219	! dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
220	! ENDDO
221	! ENDIF
222	! fin changement 10 07 96
223	!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
224
225	! calcul des pentes limitees
226
227	DO ij=iip2,ip1jm
228	IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
229	dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
230	ELSE
231	dyq(ij,l)=0.
232	ENDIF
233	ENDDO
234
235	ENDDO
236
237	DO l=1,llm
238	DO ij=1,ip1jm
239	IF(masse_adv_v(ij,l).gt.0) THEN
240	qbyv(ij,l)=q(ij+iip1,l)+dyq(ij+iip1,l)* &
241	& 0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l))
242	ELSE
243	qbyv(ij,l)=q(ij,l)-dyq(ij,l)* &
244	& 0.5*(1.+masse_adv_v(ij,l)/masse(ij,l))
245	ENDIF
246	qbyv(ij,l)=masse_adv_v(ij,l)*qbyv(ij,l)
247	ENDDO
248	ENDDO
249
250
251	DO l=1,llm
252	DO ij=iip2,ip1jm
253	newmasse=masse(ij,l) &
254	& +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
255	q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l)) &
256	& /newmasse
257	masse(ij,l)=newmasse
258	ENDDO
259	!.-. ancienne version
260	! convpn=SSUM(iim,qbyv(1,l),1)/apoln
261	! convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
262
263	convpn=SSUM(iim,qbyv(1,l),1)
264	convmpn=ssum(iim,masse_adv_v(1,l),1)
265	massepn=ssum(iim,masse(1,l),1)
266	qpn=0.
267	do ij=1,iim
268	qpn=qpn+masse(ij,l)*q(ij,l)
269	enddo
270	qpn=(qpn+convpn)/(massepn+convmpn)
271	do ij=1,iip1
272	q(ij,l)=qpn
273	enddo
274
275	! convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
276	! convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)/apols
277
278	convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
279	convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
280	masseps=ssum(iim, masse(ip1jm+1,l),1)
281	qps=0.
282	do ij = ip1jm+1,ip1jmp1-1
283	qps=qps+masse(ij,l)*q(ij,l)
284	enddo
285	qps=(qps+convps)/(masseps+convmps)
286	do ij=ip1jm+1,ip1jmp1
287	q(ij,l)=qps
288	enddo
289
290	!.-. fin ancienne version
291
292	!._. nouvelle version
293	! convpn=SSUM(iim,qbyv(1,l),1)
294	! convmpn=ssum(iim,masse_adv_v(1,l),1)
295	! oldmasse=ssum(iim,masse(1,l),1)
296	! newmasse=oldmasse+convmpn
297	! newq=(q(1,l)*oldmasse+convpn)/newmasse
298	! newmasse=newmasse/apoln
299	! DO ij = 1,iip1
300	! q(ij,l)=newq
301	! masse(ij,l)=newmasse*aire(ij)
302	! ENDDO
303	! convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
304	! convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
305	! oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
306	! newmasse=oldmasse+convmps
307	! newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
308	! newmasse=newmasse/apols
309	! DO ij = ip1jm+1,ip1jmp1
310	! q(ij,l)=newq
311	! masse(ij,l)=newmasse*aire(ij)
312	! ENDDO
313	!._. fin nouvelle version
314	ENDDO
315
316	RETURN
317	END