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	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	use disvert_m
14	use conf_gcm_m
15	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