dyn3d/Vlsplt/vlx.f

module vlx_m

  IMPLICIT NONE

contains

  SUBROUTINE vlx(q, pente_max, masse, u_m)

    ! Authors: P. Le Van, F. Hourdin, F. Forget

    ! Sch\'ema d'advection "pseudo-amont".

    use dimensions, only: iim, llm
    use paramet_m, only: ip1jmp1, iip1, iip2, ip1jm

    REAL, intent(inout):: q(ip1jmp1, llm)
    REAL, intent(in):: pente_max
    real, intent(inout):: masse(ip1jmp1, llm)
    REAL, intent(in):: u_m(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(iip2:ip1jm)
    REAL zz(ip1jmp1)
    REAL adxqu(iip2:ip1jm), dxqmax(ip1jmp1, llm)
    REAL u_mq(ip1jmp1, llm)

    !-----------------------------------------------------

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

    IF (pente_max > - 1e-5) 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)
          ENDDO
          DO ij = iip1 + iip1, ip1jm, iip1
             dxqu(ij) = dxqu(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))
          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) > 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
    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) > 0) THEN
                dxq(ij, l) = zz(ij) / (dxqu(ij - 1) + dxqu(ij))
             ELSE
                ! extremum local
                dxq(ij, l) = 0.
             ENDIF
          ENDDO
       ENDDO
    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.

    DO l = 1, llm
       DO ij = iip2, ip1jm - 1
          IF (u_m(ij, l) > 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

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

    ! 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 > 0) THEN
       DO l = 1, llm
          IF (nl(l) > 0) THEN
             iju = 0
             ! indicage des mailles concernees par le traitement special
             DO ij = iip2, ip1jm
                IF (iadvplus(ij, l) == 1.and.mod(ij, iip1) /= 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 > 0.) THEN
                   ijq = ij
                   i = ijq - (j - 1) * iip1
                   ! accumulation pour les mailles completements advectees
                   do while(zu_m > 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 > 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

    ! 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

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

  END SUBROUTINE vlx

end module vlx_m
1	guez	157	module vlx_m
2	guez	31
3	guez	157	IMPLICIT NONE
4	guez	31
5	guez	157	contains
6	guez	31
7	guez	157	SUBROUTINE vlx(q, pente_max, masse, u_m)
8	guez	31
9	guez	157	! Authors: P. Le Van, F. Hourdin, F. Forget
10	guez	31
11	guez	157	! Sch\'ema d'advection "pseudo-amont".
12	guez	31
13	guez	265	use dimensions, only: iim, llm
14	guez	157	use paramet_m, only: ip1jmp1, iip1, iip2, ip1jm
15	guez	31
16	guez	157	REAL, intent(inout):: q(ip1jmp1, llm)
17			REAL, intent(in):: pente_max
18			real, intent(inout):: masse(ip1jmp1, llm)
19			REAL, intent(in):: u_m(ip1jmp1, llm)
20	guez	31
21	guez	157	! Local:
22	guez	31
23	guez	157	INTEGER ij, l, j, i, iju, ijq, indu(ip1jmp1), niju
24			INTEGER n0, iadvplus(ip1jmp1, llm), nl(llm)
25			REAL new_m, zu_m, zdum(ip1jmp1, llm)
26			REAL dxq(ip1jmp1, llm), dxqu(iip2:ip1jm)
27			REAL zz(ip1jmp1)
28			REAL adxqu(iip2:ip1jm), dxqmax(ip1jmp1, llm)
29			REAL u_mq(ip1jmp1, llm)
30	guez	31
31	guez	157	!-----------------------------------------------------
32	guez	31
33	guez	157	! calcul de la pente a droite et a gauche de la maille
34	guez	31
35	guez	157	IF (pente_max > - 1e-5) THEN
36			! calcul des pentes avec limitation, Van Leer scheme I:
37	guez	31
38	guez	157	! calcul de la pente aux points u
39			DO l = 1, llm
40			DO ij = iip2, ip1jm - 1
41			dxqu(ij) = q(ij + 1, l) - q(ij, l)
42			ENDDO
43			DO ij = iip1 + iip1, ip1jm, iip1
44			dxqu(ij) = dxqu(ij - iim)
45			ENDDO
46	guez	31
47	guez	157	DO ij = iip2, ip1jm
48			adxqu(ij) = abs(dxqu(ij))
49			ENDDO
50	guez	31
51	guez	157	! calcul de la pente maximum dans la maille en valeur absolue
52	guez	31
53	guez	157	DO ij = iip2 + 1, ip1jm
54			dxqmax(ij, l) = pente_max * min(adxqu(ij - 1), adxqu(ij))
55			ENDDO
56	guez	31
57	guez	157	DO ij = iip1 + iip1, ip1jm, iip1
58			dxqmax(ij - iim, l) = dxqmax(ij, l)
59			ENDDO
60	guez	31
61	guez	157	DO ij = iip2 + 1, ip1jm
62			IF (dxqu(ij - 1) * dxqu(ij) > 0.) THEN
63			dxq(ij, l) = dxqu(ij - 1) + dxqu(ij)
64			ELSE
65			! extremum local
66			dxq(ij, l) = 0.
67			ENDIF
68			dxq(ij, l) = 0.5 * dxq(ij, l)
69			dxq(ij, l) = sign(min(abs(dxq(ij, l)), dxqmax(ij, l)), dxq(ij, l))
70			ENDDO
71			ENDDO
72			ELSE
73			! Pentes produits:
74	guez	31
75	guez	157	DO l = 1, llm
76			DO ij = iip2, ip1jm - 1
77			dxqu(ij) = q(ij + 1, l) - q(ij, l)
78			ENDDO
79			DO ij = iip1 + iip1, ip1jm, iip1
80			dxqu(ij) = dxqu(ij - iim)
81			ENDDO
82	guez	31
83	guez	157	DO ij = iip2 + 1, ip1jm
84			zz(ij) = dxqu(ij - 1) * dxqu(ij)
85			zz(ij) = zz(ij) + zz(ij)
86			IF (zz(ij) > 0) THEN
87			dxq(ij, l) = zz(ij) / (dxqu(ij - 1) + dxqu(ij))
88			ELSE
89			! extremum local
90			dxq(ij, l) = 0.
91			ENDIF
92			ENDDO
93			ENDDO
94			ENDIF
95	guez	31
96	guez	157	! bouclage de la pente en iip1:
97	guez	31
98	guez	157	DO l = 1, llm
99			DO ij = iip1 + iip1, ip1jm, iip1
100			dxq(ij - iim, l) = dxq(ij, l)
101			ENDDO
102			DO ij = 1, ip1jmp1
103			iadvplus(ij, l) = 0
104			ENDDO
105			ENDDO
106	guez	31
107	guez	157	! calcul des flux a gauche et a droite
108	guez	31
109	guez	157	! on cumule le flux correspondant a toutes les mailles dont la masse
110			! au travers de la paroi pENDant le pas de temps.
111	guez	31
112	guez	157	DO l = 1, llm
113			DO ij = iip2, ip1jm - 1
114			IF (u_m(ij, l) > 0.) THEN
115			zdum(ij, l) = 1. - u_m(ij, l) / masse(ij, l)
116			u_mq(ij, l) = u_m(ij, l) &
117			* (q(ij, l) + 0.5 * zdum(ij, l) * dxq(ij, l))
118	guez	31	ELSE
119	guez	157	zdum(ij, l) = 1. + u_m(ij, l) / masse(ij + 1, l)
120			u_mq(ij, l) = u_m(ij, l) &
121			* (q(ij + 1, l) - 0.5 * zdum(ij, l) * dxq(ij + 1, l))
122	guez	31	ENDIF
123			ENDDO
124	guez	157	ENDDO
125	guez	31
126	guez	157	! detection des points ou on advecte plus que la masse de la
127			! maille
128			DO l = 1, llm
129			DO ij = iip2, ip1jm - 1
130			IF (zdum(ij, l) < 0) THEN
131			iadvplus(ij, l) = 1
132			u_mq(ij, l) = 0.
133			ENDIF
134	guez	31	ENDDO
135	guez	157	ENDDO
136	guez	31
137	guez	157	DO l = 1, llm
138			DO ij = iip1 + iip1, ip1jm, iip1
139			iadvplus(ij, l) = iadvplus(ij - iim, l)
140			ENDDO
141			ENDDO
142	guez	31
143	guez	157	! traitement special pour le cas ou on advecte en longitude plus
144			! que le contenu de la maille. cette partie est mal vectorisee.
145	guez	31
146	guez	157	! calcul du nombre de maille sur lequel on advecte plus que la maille.
147	guez	31
148	guez	157	n0 = 0
149			DO l = 1, llm
150			nl(l) = 0
151			DO ij = iip2, ip1jm
152			nl(l) = nl(l) + iadvplus(ij, l)
153			ENDDO
154			n0 = n0 + nl(l)
155			ENDDO
156	guez	31
157	guez	157	IF (n0 > 0) THEN
158			DO l = 1, llm
159			IF (nl(l) > 0) THEN
160			iju = 0
161			! indicage des mailles concernees par le traitement special
162			DO ij = iip2, ip1jm
163			IF (iadvplus(ij, l) == 1.and.mod(ij, iip1) /= 0) THEN
164			iju = iju + 1
165			indu(iju) = ij
166			ENDIF
167			ENDDO
168			niju = iju
169	guez	31
170	guez	157	! traitement des mailles
171			DO iju = 1, niju
172			ij = indu(iju)
173			j = (ij - 1) / iip1 + 1
174			zu_m = u_m(ij, l)
175			u_mq(ij, l) = 0.
176			IF (zu_m > 0.) THEN
177			ijq = ij
178			i = ijq - (j - 1) * iip1
179			! accumulation pour les mailles completements advectees
180			do while(zu_m > masse(ijq, l))
181			u_mq(ij, l) = u_mq(ij, l) + q(ijq, l) * masse(ijq, l)
182			zu_m = zu_m - masse(ijq, l)
183			i = mod(i - 2 + iim, iim) + 1
184			ijq = (j - 1) * iip1 + i
185			ENDDO
186			! ajout de la maille non completement advectee
187			u_mq(ij, l) = u_mq(ij, l) + zu_m * (q(ijq, l) + 0.5 * (1. &
188			- zu_m / masse(ijq, l)) * dxq(ijq, l))
189			ELSE
190			ijq = ij + 1
191			i = ijq - (j - 1) * iip1
192			! accumulation pour les mailles completements advectees
193			do while(- zu_m > masse(ijq, l))
194			u_mq(ij, l) = u_mq(ij, l) - q(ijq, l) * masse(ijq, l)
195			zu_m = zu_m + masse(ijq, l)
196			i = mod(i, iim) + 1
197			ijq = (j - 1) * iip1 + i
198			ENDDO
199			! ajout de la maille non completement advectee
200			u_mq(ij, l) = u_mq(ij, l) + zu_m * (q(ijq, l) - 0.5 * (1. &
201			+ zu_m / masse(ijq, l)) * dxq(ijq, l))
202			ENDIF
203			ENDDO
204			ENDIF
205			ENDDO
206			ENDIF
207	guez	31
208	guez	157	! bouclage en latitude
209			DO l = 1, llm
210			DO ij = iip1 + iip1, ip1jm, iip1
211			u_mq(ij, l) = u_mq(ij - iim, l)
212			ENDDO
213			ENDDO
214	guez	31
215	guez	157	! calcul des tENDances
216	guez	31
217	guez	157	DO l = 1, llm
218			DO ij = iip2 + 1, ip1jm
219			new_m = masse(ij, l) + u_m(ij - 1, l) - u_m(ij, l)
220			q(ij, l) = (q(ij, l) * masse(ij, l) + u_mq(ij - 1, l) &
221			- u_mq(ij, l)) / new_m
222			masse(ij, l) = new_m
223			ENDDO
224	guez	31
225	guez	157	DO ij = iip1 + iip1, ip1jm, iip1
226			q(ij - iim, l) = q(ij, l)
227			masse(ij - iim, l) = masse(ij, l)
228			ENDDO
229			ENDDO
230	guez	31
231	guez	157	END SUBROUTINE vlx
232	guez	31
233	guez	157	end module vlx_m