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