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 massbar_m, only: massbar
    USE paramet_m, ONLY : iip1, iip2, ijmllm, ijp1llm, ip1jm, ip1jmp1, jjp1, &
         llmp1

    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:: dum = .TRUE., 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(massem, 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ées sur la masse
       CALL massbar(massem, massebx, masseby)

       ! Appel des sous programmes d'advection

       DO iq = 1, nqmx
          IF (iadv(iq)==0) CYCLE

          ! Schema de Van Leer I MUSCL

          IF (iadv(iq)==10) THEN
             CALL vlsplt(q(:, :, iq), 2., massem, wg, pbarug, pbarvg, dtvr)
             ! Schema "pseudo amont" + test sur humidite specifique
             ! pour la vapeur d'eau. F. Codron
          ELSE IF (iadv(iq)==14) THEN
             CALL vlspltqs(q(1, 1, 1), 2., massem, wg, pbarug, pbarvg, dtvr, &
                  p, pk, teta)
             ! Schema de Frederic Hourdin
          ELSE IF (iadv(iq)==12) THEN
             ! 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
          ELSE IF (iadv(iq)==13) THEN
             ! 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
             ! Schema de pente SLOPES
          ELSE IF (iadv(iq)==20) THEN
             CALL pentes_ini(q(1, 1, iq), wg, massem, pbarug, pbarvg, 0)
             ! Schema de Prather
          ELSE IF (iadv(iq)==30) THEN
             ! 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)
             ! Schemas PPM Lin et Rood
          ELSE IF (iadv(iq)==11 .OR. (iadv(iq)>=16 .AND. iadv(iq)<=18)) THEN
             ! 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, dum, 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, dum, 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, dum, 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, dum, 220.)
                END IF
             END DO

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

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

  END SUBROUTINE advtrac

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