trunk/dyn3d/advtrac.f

module advtrac_m

  IMPLICIT NONE

contains

  SUBROUTINE advtrac(pbaru, pbarv, p, masse, q, iapptrac, teta, pk)

    ! From dyn3d/advtrac.F, version 1.4 2005/04/13 08:58:34
    ! Author: F. Hourdin

    USE comconst, ONLY : dtvr
    USE conf_gcm_m, ONLY : iapp_tracvl
    USE dimens_m, ONLY : iim, jjm, llm, nqmx
    USE iniadvtrac_m, ONLY : iadv
    use interpre_m, only: interpre
    use massbar_m, only: massbar
    USE paramet_m, ONLY : iip1, iip2, ijmllm, ijp1llm, ip1jm, ip1jmp1, jjp1, &
         llmp1
    use ppm3d_m, only: ppm3d
    use vlsplt_m, only: vlsplt
    use vlspltqs_m, only: vlspltqs

    REAL, intent(in):: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
    REAL, intent(in):: p(ip1jmp1, llmp1)
    real, intent(in):: masse(ip1jmp1, llm)
    REAL, intent(inout):: q(ip1jmp1, llm, nqmx)
    INTEGER, intent(out):: iapptrac
    real, intent(in):: teta(ip1jmp1, llm)
    REAL, intent(in):: pk(ip1jmp1, llm)

    ! Variables locales

    REAL massebx(ip1jmp1, llm), masseby(ip1jm, llm)
    REAL, save:: pbaruc(ip1jmp1, llm), pbarvc(ip1jm, llm)
    REAL, save:: massem(ip1jmp1, llm)
    real zdp(ip1jmp1)
    REAL pbarug(ip1jmp1, llm), pbarvg(ip1jm, llm), wg(ip1jmp1, llm)

    INTEGER:: iadvtr = 0
    INTEGER ij, l, iq
    REAL zdpmin, zdpmax
    EXTERNAL minmax

    ! Rajouts pour PPM

    INTEGER indice, n
    ! Pas de temps adaptatif pour que CFL < 1 
    REAL dtbon
    REAL cflmaxz ! CFL maximum
    real aaa, bbb
    REAL psppm(iim, jjp1) ! pression au sol 
    REAL unatppm(iim, jjp1, llm), vnatppm(iim, jjp1, llm)
    REAL qppm(iim*jjp1, llm, nqmx)
    REAL fluxwppm(iim, jjp1, llm)
    REAL apppm(llmp1), bpppm(llmp1)
    LOGICAL:: fill = .TRUE.

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

    IF (iadvtr==0) THEN
       CALL initial0(ijp1llm, pbaruc)
       CALL initial0(ijmllm, pbarvc)
    END IF

    ! accumulation des flux de masse horizontaux
    DO l = 1, llm
       DO ij = 1, ip1jmp1
          pbaruc(ij, l) = pbaruc(ij, l) + pbaru(ij, l)
       END DO
       DO ij = 1, ip1jm
          pbarvc(ij, l) = pbarvc(ij, l) + pbarv(ij, l)
       END DO
    END DO

    ! selection de la masse instantannee des mailles avant le transport.
    IF (iadvtr==0) massem = masse

    iadvtr = iadvtr + 1
    iapptrac = iadvtr

    ! Test pour savoir si on advecte a ce pas de temps
    IF (iadvtr == iapp_tracvl) THEN
       ! traitement des flux de masse avant advection.
       ! 1. calcul de w
       ! 2. groupement des mailles pres du pole.

       CALL groupe(pbaruc, pbarvc, pbarug, pbarvg, wg)

       ! test sur l'eventuelle creation de valeurs negatives de la masse
       DO l = 1, llm - 1
          DO ij = iip2 + 1, ip1jm
             zdp(ij) = pbarug(ij-1, l) - pbarug(ij, l) - pbarvg(ij-iip1, l) + &
                  pbarvg(ij, l) + wg(ij, l+1) - wg(ij, l)
          END DO
          CALL scopy(jjm-1, zdp(iip1+iip1), iip1, zdp(iip2), iip1)
          DO ij = iip2, ip1jm
             zdp(ij) = zdp(ij)*dtvr/massem(ij, l)
          END DO

          CALL minmax(ip1jm-iip1, zdp(iip2), zdpmin, zdpmax)

          IF (max(abs(zdpmin), abs(zdpmax))>0.5) THEN
             PRINT *, 'WARNING DP/P l=', l, ' MIN:', zdpmin, ' MAX:', zdpmax
          END IF
       END DO

       ! Advection proprement dite

       ! Calcul des moyennes bas\'ees sur la masse
       CALL massbar(massem, massebx, masseby)

       ! Appel des sous programmes d'advection

       DO iq = 1, nqmx
          select case (iadv(iq))
          case (10)
             ! Schema de Van Leer I MUSCL
             CALL vlsplt(q(:, :, iq), 2., massem, wg, pbarug, pbarvg, dtvr)
          case (14)
             ! Schema "pseudo amont" + test sur humidite specifique
             ! pour la vapeur d'eau. F. Codron
             CALL vlspltqs(q(1, 1, 1), 2., massem, wg, pbarug, pbarvg, dtvr, &
                  p, pk, teta)
          case (12)
             ! Schema de Frederic Hourdin
             ! Pas de temps adaptatif
             CALL adaptdt(dtbon, n, pbarug, massem)
             IF (n>1) THEN
                WRITE (*, *) 'WARNING horizontal dt=', dtbon, 'dtvr=', dtvr, &
                     'n=', n
             END IF
             DO indice = 1, n
                CALL advn(q(1, 1, iq), massem, wg, pbarug, pbarvg, dtbon, 1)
             END DO
          case (13)
             ! Pas de temps adaptatif
             CALL adaptdt(dtbon, n, pbarug, massem)
             IF (n>1) THEN
                WRITE (*, *) 'WARNING horizontal dt=', dtbon, 'dtvr=', dtvr, &
                     'n=', n
             END IF
             DO indice = 1, n
                CALL advn(q(1, 1, iq), massem, wg, pbarug, pbarvg, dtbon, 2)
             END DO
          case (20)
             ! Schema de pente SLOPES
             CALL pentes_ini(q(1, 1, iq), wg, massem, pbarug, pbarvg, 0)
          case (30)
             ! Schema de Prather
             ! Pas de temps adaptatif
             CALL adaptdt(dtbon, n, pbarug, massem)
             IF (n>1) THEN
                WRITE (*, *) 'WARNING horizontal dt=', dtbon, 'dtvr=', dtvr, &
                     'n=', n
             END IF
             CALL prather(q(1, 1, iq), wg, massem, pbarug, pbarvg, n, dtbon)
          case (11, 16:18)
             ! Schemas PPM Lin et Rood
             ! Test sur le flux horizontal
             ! Pas de temps adaptatif
             CALL adaptdt(dtbon, n, pbarug, massem)
             IF (n>1) THEN
                WRITE (*, *) 'WARNING horizontal dt=', dtbon, 'dtvr=', dtvr, &
                     'n=', n
             END IF
             ! Test sur le flux vertical
             cflmaxz = 0.
             DO l = 2, llm
                DO ij = iip2, ip1jm
                   aaa = wg(ij, l)*dtvr/massem(ij, l)
                   cflmaxz = max(cflmaxz, aaa)
                   bbb = -wg(ij, l)*dtvr/massem(ij, l-1)
                   cflmaxz = max(cflmaxz, bbb)
                END DO
             END DO
             IF (cflmaxz>=1) THEN
                WRITE (*, *) 'WARNING vertical', 'CFLmaxz=', cflmaxz
             END IF

             ! Ss-prg interface LMDZ.4->PPM3d
             CALL interpre(q(1, 1, iq), qppm(1, 1, iq), wg, fluxwppm, massem, &
                  apppm, bpppm, massebx, masseby, pbarug, pbarvg, unatppm, &
                  vnatppm, psppm)

             DO indice = 1, n
                ! VL (version PPM) horiz. et PPM vert.
                IF (iadv(iq)==11) THEN
                   ! Ss-prg PPM3d de Lin
                   CALL ppm3d(1, qppm(1, 1, iq), psppm, psppm, unatppm, &
                        vnatppm, fluxwppm, dtbon, 2, 2, 2, 1, iim, jjp1, 2, &
                        llm, apppm, bpppm, 0.01, 6400000, fill, 220.)
                   ! Monotonic PPM
                ELSE IF (iadv(iq)==16) THEN
                   ! Ss-prg PPM3d de Lin
                   CALL ppm3d(1, qppm(1, 1, iq), psppm, psppm, unatppm, &
                        vnatppm, fluxwppm, dtbon, 3, 3, 3, 1, iim, jjp1, 2, &
                        llm, apppm, bpppm, 0.01, 6400000, fill, 220.)
                   ! Semi Monotonic PPM
                ELSE IF (iadv(iq)==17) THEN
                   ! Ss-prg PPM3d de Lin
                   CALL ppm3d(1, qppm(1, 1, iq), psppm, psppm, unatppm, &
                        vnatppm, fluxwppm, dtbon, 4, 4, 4, 1, iim, jjp1, 2, &
                        llm, apppm, bpppm, 0.01, 6400000, fill, 220.)
                   ! Positive Definite PPM
                ELSE IF (iadv(iq)==18) THEN
                   ! Ss-prg PPM3d de Lin
                   CALL ppm3d(1, qppm(1, 1, iq), psppm, psppm, unatppm, &
                        vnatppm, fluxwppm, dtbon, 5, 5, 5, 1, iim, jjp1, 2, &
                        llm, apppm, bpppm, 0.01, 6400000, fill, 220.)
                END IF
             END DO

             ! Ss-prg interface PPM3d-LMDZ.4
             CALL interpost(q(1, 1, iq), qppm(1, 1, iq))
          END select
       END DO

       ! on reinitialise a zero les flux de masse cumules
       iadvtr = 0
    END IF

  END SUBROUTINE advtrac

end module advtrac_m
1	guez	40	module advtrac_m
2	guez	3
3	guez	40	IMPLICIT NONE
4	guez	3
5	guez	40	contains
6	guez	23
7	guez	40	SUBROUTINE advtrac(pbaru, pbarv, p, masse, q, iapptrac, teta, pk)
8	guez	3
9	guez	40	! From dyn3d/advtrac.F, version 1.4 2005/04/13 08:58:34
10			! Author: F. Hourdin
11	guez	3
12	guez	71	USE comconst, ONLY : dtvr
13			USE conf_gcm_m, ONLY : iapp_tracvl
14	guez	40	USE dimens_m, ONLY : iim, jjm, llm, nqmx
15	guez	71	USE iniadvtrac_m, ONLY : iadv
16	guez	108	use interpre_m, only: interpre
17	guez	91	use massbar_m, only: massbar
18	guez	40	USE paramet_m, ONLY : iip1, iip2, ijmllm, ijp1llm, ip1jm, ip1jmp1, jjp1, &
19			llmp1
20	guez	167	use ppm3d_m, only: ppm3d
21	guez	157	use vlsplt_m, only: vlsplt
22	guez	108	use vlspltqs_m, only: vlspltqs
23	guez	31
24	guez	71	REAL, intent(in):: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
25			REAL, intent(in):: p(ip1jmp1, llmp1)
26			real, intent(in):: masse(ip1jmp1, llm)
27	guez	40	REAL, intent(inout):: q(ip1jmp1, llm, nqmx)
28	guez	71	INTEGER, intent(out):: iapptrac
29	guez	44	real, intent(in):: teta(ip1jmp1, llm)
30	guez	71	REAL, intent(in):: pk(ip1jmp1, llm)
31	guez	3
32	guez	40	! Variables locales
33	guez	3
34	guez	40	REAL massebx(ip1jmp1, llm), masseby(ip1jm, llm)
35			REAL, save:: pbaruc(ip1jmp1, llm), pbarvc(ip1jm, llm)
36			REAL, save:: massem(ip1jmp1, llm)
37			real zdp(ip1jmp1)
38			REAL pbarug(ip1jmp1, llm), pbarvg(ip1jm, llm), wg(ip1jmp1, llm)
39	guez	31
40	guez	40	INTEGER:: iadvtr = 0
41			INTEGER ij, l, iq
42			REAL zdpmin, zdpmax
43			EXTERNAL minmax
44	guez	31
45	guez	40	! Rajouts pour PPM
46	guez	3
47	guez	40	INTEGER indice, n
48			! Pas de temps adaptatif pour que CFL < 1
49			REAL dtbon
50			REAL cflmaxz ! CFL maximum
51			real aaa, bbb
52			REAL psppm(iim, jjp1) ! pression au sol
53			REAL unatppm(iim, jjp1, llm), vnatppm(iim, jjp1, llm)
54			REAL qppm(iim*jjp1, llm, nqmx)
55			REAL fluxwppm(iim, jjp1, llm)
56			REAL apppm(llmp1), bpppm(llmp1)
57	guez	168	LOGICAL:: fill = .TRUE.
58	guez	3
59	guez	40	!-----------------------------------------------------------
60	guez	3
61	guez	40	IF (iadvtr==0) THEN
62			CALL initial0(ijp1llm, pbaruc)
63			CALL initial0(ijmllm, pbarvc)
64			END IF
65	guez	3
66	guez	40	! accumulation des flux de masse horizontaux
67			DO l = 1, llm
68			DO ij = 1, ip1jmp1
69			pbaruc(ij, l) = pbaruc(ij, l) + pbaru(ij, l)
70			END DO
71			DO ij = 1, ip1jm
72			pbarvc(ij, l) = pbarvc(ij, l) + pbarv(ij, l)
73			END DO
74			END DO
75	guez	3
76	guez	40	! selection de la masse instantannee des mailles avant le transport.
77	guez	71	IF (iadvtr==0) massem = masse
78	guez	3
79	guez	40	iadvtr = iadvtr + 1
80			iapptrac = iadvtr
81	guez	3
82	guez	40	! Test pour savoir si on advecte a ce pas de temps
83	guez	71	IF (iadvtr == iapp_tracvl) THEN
84	guez	40	! traitement des flux de masse avant advection.
85			! 1. calcul de w
86			! 2. groupement des mailles pres du pole.
87	guez	3
88	guez	150	CALL groupe(pbaruc, pbarvc, pbarug, pbarvg, wg)
89	guez	3
90	guez	40	! test sur l'eventuelle creation de valeurs negatives de la masse
91			DO l = 1, llm - 1
92			DO ij = iip2 + 1, ip1jm
93			zdp(ij) = pbarug(ij-1, l) - pbarug(ij, l) - pbarvg(ij-iip1, l) + &
94			pbarvg(ij, l) + wg(ij, l+1) - wg(ij, l)
95			END DO
96			CALL scopy(jjm-1, zdp(iip1+iip1), iip1, zdp(iip2), iip1)
97			DO ij = iip2, ip1jm
98			zdp(ij) = zdp(ij)*dtvr/massem(ij, l)
99			END DO
100	guez	3
101	guez	40	CALL minmax(ip1jm-iip1, zdp(iip2), zdpmin, zdpmax)
102	guez	3
103	guez	40	IF (max(abs(zdpmin), abs(zdpmax))>0.5) THEN
104			PRINT *, 'WARNING DP/P l=', l, ' MIN:', zdpmin, ' MAX:', zdpmax
105			END IF
106			END DO
107	guez	3
108	guez	40	! Advection proprement dite
109	guez	3
110	guez	150	! Calcul des moyennes bas\'ees sur la masse
111	guez	40	CALL massbar(massem, massebx, masseby)
112	guez	31
113	guez	40	! Appel des sous programmes d'advection
114	guez	31
115	guez	40	DO iq = 1, nqmx
116	guez	157	select case (iadv(iq))
117			case (10)
118			! Schema de Van Leer I MUSCL
119	guez	40	CALL vlsplt(q(:, :, iq), 2., massem, wg, pbarug, pbarvg, dtvr)
120	guez	157	case (14)
121	guez	40	! Schema "pseudo amont" + test sur humidite specifique
122			! pour la vapeur d'eau. F. Codron
123			CALL vlspltqs(q(1, 1, 1), 2., massem, wg, pbarug, pbarvg, dtvr, &
124			p, pk, teta)
125	guez	157	case (12)
126	guez	40	! Schema de Frederic Hourdin
127			! Pas de temps adaptatif
128	guez	104	CALL adaptdt(dtbon, n, pbarug, massem)
129	guez	40	IF (n>1) THEN
130			WRITE (, ) 'WARNING horizontal dt=', dtbon, 'dtvr=', dtvr, &
131			'n=', n
132			END IF
133			DO indice = 1, n
134			CALL advn(q(1, 1, iq), massem, wg, pbarug, pbarvg, dtbon, 1)
135			END DO
136	guez	157	case (13)
137	guez	40	! Pas de temps adaptatif
138	guez	104	CALL adaptdt(dtbon, n, pbarug, massem)
139	guez	40	IF (n>1) THEN
140			WRITE (, ) 'WARNING horizontal dt=', dtbon, 'dtvr=', dtvr, &
141			'n=', n
142			END IF
143			DO indice = 1, n
144			CALL advn(q(1, 1, iq), massem, wg, pbarug, pbarvg, dtbon, 2)
145			END DO
146	guez	157	case (20)
147	guez	40	! Schema de pente SLOPES
148			CALL pentes_ini(q(1, 1, iq), wg, massem, pbarug, pbarvg, 0)
149	guez	157	case (30)
150	guez	40	! Schema de Prather
151			! Pas de temps adaptatif
152	guez	104	CALL adaptdt(dtbon, n, pbarug, massem)
153	guez	40	IF (n>1) THEN
154			WRITE (, ) 'WARNING horizontal dt=', dtbon, 'dtvr=', dtvr, &
155			'n=', n
156			END IF
157			CALL prather(q(1, 1, iq), wg, massem, pbarug, pbarvg, n, dtbon)
158	guez	157	case (11, 16:18)
159	guez	40	! Schemas PPM Lin et Rood
160			! Test sur le flux horizontal
161			! Pas de temps adaptatif
162	guez	104	CALL adaptdt(dtbon, n, pbarug, massem)
163	guez	40	IF (n>1) THEN
164			WRITE (, ) 'WARNING horizontal dt=', dtbon, 'dtvr=', dtvr, &
165			'n=', n
166			END IF
167			! Test sur le flux vertical
168			cflmaxz = 0.
169			DO l = 2, llm
170			DO ij = iip2, ip1jm
171			aaa = wg(ij, l)*dtvr/massem(ij, l)
172			cflmaxz = max(cflmaxz, aaa)
173			bbb = -wg(ij, l)*dtvr/massem(ij, l-1)
174			cflmaxz = max(cflmaxz, bbb)
175			END DO
176			END DO
177			IF (cflmaxz>=1) THEN
178			WRITE (, ) 'WARNING vertical', 'CFLmaxz=', cflmaxz
179			END IF
180	guez	3
181	guez	40	! Ss-prg interface LMDZ.4->PPM3d
182			CALL interpre(q(1, 1, iq), qppm(1, 1, iq), wg, fluxwppm, massem, &
183			apppm, bpppm, massebx, masseby, pbarug, pbarvg, unatppm, &
184			vnatppm, psppm)
185	guez	31
186	guez	40	DO indice = 1, n
187			! VL (version PPM) horiz. et PPM vert.
188			IF (iadv(iq)==11) THEN
189			! Ss-prg PPM3d de Lin
190			CALL ppm3d(1, qppm(1, 1, iq), psppm, psppm, unatppm, &
191			vnatppm, fluxwppm, dtbon, 2, 2, 2, 1, iim, jjp1, 2, &
192	guez	167	llm, apppm, bpppm, 0.01, 6400000, fill, 220.)
193	guez	40	! Monotonic PPM
194			ELSE IF (iadv(iq)==16) THEN
195			! Ss-prg PPM3d de Lin
196			CALL ppm3d(1, qppm(1, 1, iq), psppm, psppm, unatppm, &
197			vnatppm, fluxwppm, dtbon, 3, 3, 3, 1, iim, jjp1, 2, &
198	guez	167	llm, apppm, bpppm, 0.01, 6400000, fill, 220.)
199	guez	40	! Semi Monotonic PPM
200			ELSE IF (iadv(iq)==17) THEN
201			! Ss-prg PPM3d de Lin
202			CALL ppm3d(1, qppm(1, 1, iq), psppm, psppm, unatppm, &
203			vnatppm, fluxwppm, dtbon, 4, 4, 4, 1, iim, jjp1, 2, &
204	guez	167	llm, apppm, bpppm, 0.01, 6400000, fill, 220.)
205	guez	40	! Positive Definite PPM
206			ELSE IF (iadv(iq)==18) THEN
207			! Ss-prg PPM3d de Lin
208			CALL ppm3d(1, qppm(1, 1, iq), psppm, psppm, unatppm, &
209			vnatppm, fluxwppm, dtbon, 5, 5, 5, 1, iim, jjp1, 2, &
210	guez	167	llm, apppm, bpppm, 0.01, 6400000, fill, 220.)
211	guez	40	END IF
212			END DO
213	guez	3
214	guez	40	! Ss-prg interface PPM3d-LMDZ.4
215			CALL interpost(q(1, 1, iq), qppm(1, 1, iq))
216	guez	157	END select
217	guez	40	END DO
218	guez	3
219	guez	40	! on reinitialise a zero les flux de masse cumules
220			iadvtr = 0
221			END IF
222
223			END SUBROUTINE advtrac
224
225			end module advtrac_m