trunk/dyn3d/vlxqs.f

      SUBROUTINE vlxqs(q,pente_max,masse,u_m,qsat)
!
!     Auteurs:   P.Le Van, F.Hourdin, F.Forget
!
!    ********************************************************************
!     Shema  d'advection " pseudo amont " .
!    ********************************************************************
!
!   --------------------------------------------------------------------
      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 )
      REAL q(ip1jmp1,llm)
      REAL qsat(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 dxq(ip1jmp1,llm),dxqu(ip1jmp1)
      REAL zz(ip1jmp1)
      REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
      REAL u_mq(ip1jmp1,llm)

      Logical first,testcpu
      SAVE first,testcpu

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


      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

         ENDDO ! l=1,llm

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

!   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

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

!   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.
!   le rapport de melange de l'air advecte est min(q_vanleer, Qsat_downwind)
      DO l=1,llm
       DO ij=iip2,ip1jm-1
          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)* &
               min(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l),qsat(ij+1,l))
          ELSE
             zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
             u_mq(ij,l)=u_m(ij,l)* &
               min(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l),qsat(ij,l))
          ENDIF
       ENDDO
      ENDDO


!   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
      DO l=1,llm
       DO ij=iip1+iip1,ip1jm,iip1
          iadvplus(ij,l)=iadvplus(ij-iim,l)
       ENDDO
      ENDDO


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

!   pas d'influence de la pression saturante (pour l'instant)

!  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
         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

!  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
      ENDIF  ! n0.gt.0


!   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

!     CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
!     CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)


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