dyn3d/Vlsplt/vlx.f

      SUBROUTINE vlx(q,pente_max,masse,u_m)

!     Auteurs:   P.Le Van, F.Hourdin, F.Forget
!
!   *******************************************************
!     Shema  d'advection " pseudo amont " .
!    ************************************************************
!     nq,iq,q,pbaru,pbarv,w sont des arguments d'entree  pour le
!     s-pg ....
!
!
!   --------------------------------------------------------------
      use dimens_m
      use paramet_m
      use comconst
      use disvert_m
      use conf_gcm_m
      IMPLICIT NONE
!
!
!
!   Arguments:
!   ----------
      REAL masse(ip1jmp1,llm),pente_max
      REAL u_m( ip1jmp1,llm ),pbarv( iip1,jjm,llm)
      REAL q(ip1jmp1,llm)
      REAL w(ip1jmp1,llm)
!
!      Local
!   ---------
!
      INTEGER ij,l,j,i,iju,ijq,indu(ip1jmp1),niju
      INTEGER n0,iadvplus(ip1jmp1,llm),nl(llm)
!
      REAL new_m,zu_m,zdum(ip1jmp1,llm)
      REAL sigu(ip1jmp1),dxq(ip1jmp1,llm),dxqu(ip1jmp1)
      REAL zz(ip1jmp1)
      REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
      REAL u_mq(ip1jmp1,llm)

      Logical extremum,first,testcpu
      SAVE first,testcpu

      REAL      SSUM
      REAL temps0,temps1,temps2,temps3,temps4,temps5,second
      SAVE temps0,temps1,temps2,temps3,temps4,temps5

      REAL z1,z2,z3

      DATA first,testcpu/.true.,.false./

      IF(first) THEN
         temps1=0.
         temps2=0.
         temps3=0.
         temps4=0.
         temps5=0.
         first=.false.
      ENDIF

!   calcul de la pente a droite et a gauche de la maille


      IF (pente_max.gt.-1.e-5) THEN
!       IF (pente_max.gt.10) THEN

!   calcul des pentes avec limitation, Van Leer scheme I:
!   -----------------------------------------------------

!   calcul de la pente aux points u
         DO l = 1, llm
            DO ij=iip2,ip1jm-1
               dxqu(ij)=q(ij+1,l)-q(ij,l)
!              IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
!              sigu(ij)=u_m(ij,l)/masse(ij,l)
            ENDDO
            DO ij=iip1+iip1,ip1jm,iip1
               dxqu(ij)=dxqu(ij-iim)
!              sigu(ij)=sigu(ij-iim)
            ENDDO

            DO ij=iip2,ip1jm
               adxqu(ij)=abs(dxqu(ij))
            ENDDO

!   calcul de la pente maximum dans la maille en valeur absolue

            DO ij=iip2+1,ip1jm
               dxqmax(ij,l)=pente_max*                                  &
     &      min(adxqu(ij-1),adxqu(ij))
! limitation subtile
!    ,      min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))


            ENDDO

            DO ij=iip1+iip1,ip1jm,iip1
               dxqmax(ij-iim,l)=dxqmax(ij,l)
            ENDDO

            DO ij=iip2+1,ip1jm
               IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
                  dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
               ELSE
!   extremum local
                  dxq(ij,l)=0.
               ENDIF
               dxq(ij,l)=0.5*dxq(ij,l)
               dxq(ij,l)=                                               &
     &         sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
            ENDDO

               ! l=1,llm
         ENDDO
           ! (pente_max.lt.-1.e-5)
      ELSE

!   Pentes produits:
!   ----------------

         DO l = 1, llm
            DO ij=iip2,ip1jm-1
               dxqu(ij)=q(ij+1,l)-q(ij,l)
            ENDDO
            DO ij=iip1+iip1,ip1jm,iip1
               dxqu(ij)=dxqu(ij-iim)
            ENDDO

            DO ij=iip2+1,ip1jm
               zz(ij)=dxqu(ij-1)*dxqu(ij)
               zz(ij)=zz(ij)+zz(ij)
               IF(zz(ij).gt.0) THEN
                  dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
               ELSE
!   extremum local
                  dxq(ij,l)=0.
               ENDIF
            ENDDO

         ENDDO

            ! (pente_max.lt.-1.e-5)
      ENDIF

!   bouclage de la pente en iip1:
!   -----------------------------

      DO l=1,llm
         DO ij=iip1+iip1,ip1jm,iip1
            dxq(ij-iim,l)=dxq(ij,l)
         ENDDO
         DO ij=1,ip1jmp1
            iadvplus(ij,l)=0
         ENDDO

      ENDDO

!   calcul des flux a gauche et a droite

!   on cumule le flux correspondant a toutes les mailles dont la masse
!   au travers de la paroi pENDant le pas de temps.
!print*,'Cumule ....'

      DO l=1,llm
       DO ij=iip2,ip1jm-1
!      print*,'masse(',ij,')=',masse(ij,l)
          IF (u_m(ij,l).gt.0.) THEN
             zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
             u_mq(ij,l)=u_m(ij,l)*(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l))
          ELSE
             zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
             u_mq(ij,l)=u_m(ij,l)                                       &
     &            *(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l))
          ENDIF
       ENDDO
      ENDDO
!      stop

!      go to 9999
!   detection des points ou on advecte plus que la masse de la
!   maille
      DO l=1,llm
         DO ij=iip2,ip1jm-1
            IF(zdum(ij,l).lt.0) THEN
               iadvplus(ij,l)=1
               u_mq(ij,l)=0.
            ENDIF
         ENDDO
      ENDDO
!print*,'Ok test 1'
      DO l=1,llm
       DO ij=iip1+iip1,ip1jm,iip1
          iadvplus(ij,l)=iadvplus(ij-iim,l)
       ENDDO
      ENDDO
! print*,'Ok test 2'


!   traitement special pour le cas ou on advecte en longitude plus
!     que le
!   contenu de la maille.
!   cette partie est mal vectorisee.

!  calcul du nombre de maille sur lequel on advecte plus que la maille.

      n0=0
      DO l=1,llm
         nl(l)=0
         DO ij=iip2,ip1jm
            nl(l)=nl(l)+iadvplus(ij,l)
         ENDDO
         n0=n0+nl(l)
      ENDDO

      IF(n0.gt.0) THEN
!C      PRINT*,'Nombre de points pour lesquels on advect plus que le'
!C     &       ,'contenu de la maille : ',n0

         DO l=1,llm
            IF(nl(l).gt.0) THEN
               iju=0
!   indicage des mailles concernees par le traitement special
               DO ij=iip2,ip1jm
                  IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN
                     iju=iju+1
                     indu(iju)=ij
                  ENDIF
               ENDDO
               niju=iju
!              PRINT*,'niju,nl',niju,nl(l)

!  traitement des mailles
               DO iju=1,niju
                  ij=indu(iju)
                  j=(ij-1)/iip1+1
                  zu_m=u_m(ij,l)
                  u_mq(ij,l)=0.
                  IF(zu_m.gt.0.) THEN
                     ijq=ij
                     i=ijq-(j-1)*iip1
!   accumulation pour les mailles completements advectees
                     do while(zu_m.gt.masse(ijq,l))
                        u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
                        zu_m=zu_m-masse(ijq,l)
                        i=mod(i-2+iim,iim)+1
                        ijq=(j-1)*iip1+i
                     ENDDO
!   ajout de la maille non completement advectee
                     u_mq(ij,l)=u_mq(ij,l)+zu_m*                        &
     &               (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l))
                  ELSE
                     ijq=ij+1
                     i=ijq-(j-1)*iip1
!   accumulation pour les mailles completements advectees
                     do while(-zu_m.gt.masse(ijq,l))
                        u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
                        zu_m=zu_m+masse(ijq,l)
                        i=mod(i,iim)+1
                        ijq=(j-1)*iip1+i
                     ENDDO
!   ajout de la maille non completement advectee
                     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-              &
     &               0.5*(1.+zu_m/masse(ijq,l))*dxq(ijq,l))
                  ENDIF
               ENDDO
            ENDIF
         ENDDO
             ! n0.gt.0
      ENDIF
 9999   continue


!   bouclage en latitude
!print*,'cvant bouclage en latitude'
      DO l=1,llm
        DO ij=iip1+iip1,ip1jm,iip1
           u_mq(ij,l)=u_mq(ij-iim,l)
        ENDDO
      ENDDO


!   calcul des tENDances

      DO l=1,llm
         DO ij=iip2+1,ip1jm
            new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
            q(ij,l)=(q(ij,l)*masse(ij,l)+                               &
     &      u_mq(ij-1,l)-u_mq(ij,l))                                    &
     &      /new_m
            masse(ij,l)=new_m
         ENDDO
!   ModIF Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
         DO ij=iip1+iip1,ip1jm,iip1
            q(ij-iim,l)=q(ij,l)
            masse(ij-iim,l)=masse(ij,l)
         ENDDO
      ENDDO


      RETURN
      END
1	guez	31	SUBROUTINE vlx(q,pente_max,masse,u_m)
2
3			! Auteurs: P.Le Van, F.Hourdin, F.Forget
4			!
5			! *******************************************************
6			! Shema d'advection " pseudo amont " .
7			! ************************************************************
8			! nq,iq,q,pbaru,pbarv,w sont des arguments d'entree pour le
9			! s-pg ....
10			!
11			!
12			! --------------------------------------------------------------
13			use dimens_m
14			use paramet_m
15			use comconst
16	guez	66	use disvert_m
17	guez	57	use conf_gcm_m
18	guez	31	IMPLICIT NONE
19			!
20			!
21			!
22			! Arguments:
23			! ----------
24			REAL masse(ip1jmp1,llm),pente_max
25			REAL u_m( ip1jmp1,llm ),pbarv( iip1,jjm,llm)
26			REAL q(ip1jmp1,llm)
27			REAL w(ip1jmp1,llm)
28			!
29			! Local
30			! ---------
31			!
32			INTEGER ij,l,j,i,iju,ijq,indu(ip1jmp1),niju
33			INTEGER n0,iadvplus(ip1jmp1,llm),nl(llm)
34			!
35			REAL new_m,zu_m,zdum(ip1jmp1,llm)
36			REAL sigu(ip1jmp1),dxq(ip1jmp1,llm),dxqu(ip1jmp1)
37			REAL zz(ip1jmp1)
38			REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
39			REAL u_mq(ip1jmp1,llm)
40
41			Logical extremum,first,testcpu
42			SAVE first,testcpu
43
44			REAL SSUM
45			REAL temps0,temps1,temps2,temps3,temps4,temps5,second
46			SAVE temps0,temps1,temps2,temps3,temps4,temps5
47
48			REAL z1,z2,z3
49
50			DATA first,testcpu/.true.,.false./
51
52			IF(first) THEN
53			temps1=0.
54			temps2=0.
55			temps3=0.
56			temps4=0.
57			temps5=0.
58			first=.false.
59			ENDIF
60
61			! calcul de la pente a droite et a gauche de la maille
62
63
64			IF (pente_max.gt.-1.e-5) THEN
65			! IF (pente_max.gt.10) THEN
66
67			! calcul des pentes avec limitation, Van Leer scheme I:
68			! -----------------------------------------------------
69
70			! calcul de la pente aux points u
71			DO l = 1, llm
72			DO ij=iip2,ip1jm-1
73			dxqu(ij)=q(ij+1,l)-q(ij,l)
74			! IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
75			! sigu(ij)=u_m(ij,l)/masse(ij,l)
76			ENDDO
77			DO ij=iip1+iip1,ip1jm,iip1
78			dxqu(ij)=dxqu(ij-iim)
79			! sigu(ij)=sigu(ij-iim)
80			ENDDO
81
82			DO ij=iip2,ip1jm
83			adxqu(ij)=abs(dxqu(ij))
84			ENDDO
85
86			! calcul de la pente maximum dans la maille en valeur absolue
87
88			DO ij=iip2+1,ip1jm
89			dxqmax(ij,l)=pente_max* &
90			& min(adxqu(ij-1),adxqu(ij))
91			! limitation subtile
92			! , min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
93
94
95			ENDDO
96
97			DO ij=iip1+iip1,ip1jm,iip1
98			dxqmax(ij-iim,l)=dxqmax(ij,l)
99			ENDDO
100
101			DO ij=iip2+1,ip1jm
102			IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
103			dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
104			ELSE
105			! extremum local
106			dxq(ij,l)=0.
107			ENDIF
108			dxq(ij,l)=0.5*dxq(ij,l)
109			dxq(ij,l)= &
110			& sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
111			ENDDO
112
113			! l=1,llm
114			ENDDO
115			! (pente_max.lt.-1.e-5)
116			ELSE
117
118			! Pentes produits:
119			! ----------------
120
121			DO l = 1, llm
122			DO ij=iip2,ip1jm-1
123			dxqu(ij)=q(ij+1,l)-q(ij,l)
124			ENDDO
125			DO ij=iip1+iip1,ip1jm,iip1
126			dxqu(ij)=dxqu(ij-iim)
127			ENDDO
128
129			DO ij=iip2+1,ip1jm
130			zz(ij)=dxqu(ij-1)*dxqu(ij)
131			zz(ij)=zz(ij)+zz(ij)
132			IF(zz(ij).gt.0) THEN
133			dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
134			ELSE
135			! extremum local
136			dxq(ij,l)=0.
137			ENDIF
138			ENDDO
139
140			ENDDO
141
142			! (pente_max.lt.-1.e-5)
143			ENDIF
144
145			! bouclage de la pente en iip1:
146			! -----------------------------
147
148			DO l=1,llm
149			DO ij=iip1+iip1,ip1jm,iip1
150			dxq(ij-iim,l)=dxq(ij,l)
151			ENDDO
152			DO ij=1,ip1jmp1
153			iadvplus(ij,l)=0
154			ENDDO
155
156			ENDDO
157
158			! calcul des flux a gauche et a droite
159
160			! on cumule le flux correspondant a toutes les mailles dont la masse
161			! au travers de la paroi pENDant le pas de temps.
162			!print*,'Cumule ....'
163
164			DO l=1,llm
165			DO ij=iip2,ip1jm-1
166			! print*,'masse(',ij,')=',masse(ij,l)
167			IF (u_m(ij,l).gt.0.) THEN
168			zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
169			u_mq(ij,l)=u_m(ij,l)(q(ij,l)+0.5zdum(ij,l)*dxq(ij,l))
170			ELSE
171			zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
172			u_mq(ij,l)=u_m(ij,l) &
173			& (q(ij+1,l)-0.5zdum(ij,l)*dxq(ij+1,l))
174			ENDIF
175			ENDDO
176			ENDDO
177			! stop
178
179			! go to 9999
180			! detection des points ou on advecte plus que la masse de la
181			! maille
182			DO l=1,llm
183			DO ij=iip2,ip1jm-1
184			IF(zdum(ij,l).lt.0) THEN
185			iadvplus(ij,l)=1
186			u_mq(ij,l)=0.
187			ENDIF
188			ENDDO
189			ENDDO
190			!print*,'Ok test 1'
191			DO l=1,llm
192			DO ij=iip1+iip1,ip1jm,iip1
193			iadvplus(ij,l)=iadvplus(ij-iim,l)
194			ENDDO
195			ENDDO
196			! print*,'Ok test 2'
197
198
199			! traitement special pour le cas ou on advecte en longitude plus
200			! que le
201			! contenu de la maille.
202			! cette partie est mal vectorisee.
203
204			! calcul du nombre de maille sur lequel on advecte plus que la maille.
205
206			n0=0
207			DO l=1,llm
208			nl(l)=0
209			DO ij=iip2,ip1jm
210			nl(l)=nl(l)+iadvplus(ij,l)
211			ENDDO
212			n0=n0+nl(l)
213			ENDDO
214
215			IF(n0.gt.0) THEN
216			!C PRINT*,'Nombre de points pour lesquels on advect plus que le'
217			!C & ,'contenu de la maille : ',n0
218
219			DO l=1,llm
220			IF(nl(l).gt.0) THEN
221			iju=0
222			! indicage des mailles concernees par le traitement special
223			DO ij=iip2,ip1jm
224			IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN
225			iju=iju+1
226			indu(iju)=ij
227			ENDIF
228			ENDDO
229			niju=iju
230			! PRINT*,'niju,nl',niju,nl(l)
231
232			! traitement des mailles
233			DO iju=1,niju
234			ij=indu(iju)
235			j=(ij-1)/iip1+1
236			zu_m=u_m(ij,l)
237			u_mq(ij,l)=0.
238			IF(zu_m.gt.0.) THEN
239			ijq=ij
240			i=ijq-(j-1)*iip1
241			! accumulation pour les mailles completements advectees
242			do while(zu_m.gt.masse(ijq,l))
243			u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
244			zu_m=zu_m-masse(ijq,l)
245			i=mod(i-2+iim,iim)+1
246			ijq=(j-1)*iip1+i
247			ENDDO
248			! ajout de la maille non completement advectee
249			u_mq(ij,l)=u_mq(ij,l)+zu_m* &
250			& (q(ijq,l)+0.5(1.-zu_m/masse(ijq,l))dxq(ijq,l))
251			ELSE
252			ijq=ij+1
253			i=ijq-(j-1)*iip1
254			! accumulation pour les mailles completements advectees
255			do while(-zu_m.gt.masse(ijq,l))
256			u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
257			zu_m=zu_m+masse(ijq,l)
258			i=mod(i,iim)+1
259			ijq=(j-1)*iip1+i
260			ENDDO
261			! ajout de la maille non completement advectee
262			u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)- &
263			& 0.5(1.+zu_m/masse(ijq,l))dxq(ijq,l))
264			ENDIF
265			ENDDO
266			ENDIF
267			ENDDO
268			! n0.gt.0
269			ENDIF
270			9999 continue
271
272
273			! bouclage en latitude
274			!print*,'cvant bouclage en latitude'
275			DO l=1,llm
276			DO ij=iip1+iip1,ip1jm,iip1
277			u_mq(ij,l)=u_mq(ij-iim,l)
278			ENDDO
279			ENDDO
280
281
282			! calcul des tENDances
283
284			DO l=1,llm
285			DO ij=iip2+1,ip1jm
286			new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
287			q(ij,l)=(q(ij,l)*masse(ij,l)+ &
288			& u_mq(ij-1,l)-u_mq(ij,l)) &
289			& /new_m
290			masse(ij,l)=new_m
291			ENDDO
292			! ModIF Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
293			DO ij=iip1+iip1,ip1jm,iip1
294			q(ij-iim,l)=q(ij,l)
295			masse(ij-iim,l)=masse(ij,l)
296			ENDDO
297			ENDDO
298
299
300			RETURN
301			END