Sources/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	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 temps1,temps2,temps3,temps4,temps5
42	SAVE temps1,temps2,temps3,temps4,temps5
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