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	guez	31	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