trunk/dyn3d/advtrac.f

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 dimens_m, ONLY : iim, jjm, llm, nqmx
  USE paramet_m, ONLY : iip1, iip2, ijmllm, ijp1llm, ip1jm, ip1jmp1, jjp1, &
       llmp1
  USE comconst, ONLY : dtvr
  USE conf_gcm_m, ONLY : iapp_tracvl
  USE iniadvtrac_m, ONLY : iadv

  IMPLICIT NONE

  !     Arguments

  REAL massebx(ip1jmp1, llm), masseby(ip1jm, llm)

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

  !     Variables locales

  REAL pbaruc(ip1jmp1, llm), pbarvc(ip1jm, llm)
  REAL massem(ip1jmp1, llm), zdp(ip1jmp1)
  REAL pbarug(ip1jmp1, llm), pbarvg(ip1jm, llm), wg(ip1jmp1, llm)
  REAL cpuadv(nqmx)
  COMMON /cpuadv/cpuadv

  INTEGER iadvtr
  INTEGER ij, l, iq
  REAL zdpmin, zdpmax
  EXTERNAL minmax
  SAVE iadvtr, massem, pbaruc, pbarvc
  DATA iadvtr/0/

  !     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, fill
  DATA fill/ .TRUE./
  DATA dum/ .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) THEN
     CALL scopy(ip1jmp1*llm, masse, 1, massem, 1)
  END IF

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