Sources/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	guez	31	SUBROUTINE advtrac(pbaru, pbarv, p, masse, q, iapptrac, teta, pk)
2	guez	3
3	guez	31	! From dyn3d/advtrac.F, version 1.4 2005/04/13 08:58:34
4			! Author: F. Hourdin
5	guez	3
6	guez	23	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	guez	3	IMPLICIT NONE
14
15	guez	31	! Arguments
16	guez	3
17	guez	31	REAL massebx(ip1jmp1, llm), masseby(ip1jm, llm)
18
19	guez	3	INTEGER iapptrac
20	guez	31	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	guez	3
26			! Variables locales
27
28	guez	31	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	guez	3	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	guez	31
41	guez	3	! Rajouts pour PPM
42	guez	31
43	guez	3	INTEGER indice, n
44			! Pas de temps adaptatif pour que CFL<1
45			REAL dtbon
46	guez	31	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	guez	3	REAL apppm(llmp1), bpppm(llmp1)
53			LOGICAL dum, fill
54			DATA fill/ .TRUE./
55			DATA dum/ .TRUE./
56
57	guez	31	!-----------------------------------------------------------
58	guez	3
59			IF (iadvtr==0) THEN
60	guez	31	CALL initial0(ijp1llm, pbaruc)
61			CALL initial0(ijmllm, pbarvc)
62	guez	3	END IF
63
64			! accumulation des flux de masse horizontaux
65			DO l = 1, llm
66			DO ij = 1, ip1jmp1
67	guez	31	pbaruc(ij, l) = pbaruc(ij, l) + pbaru(ij, l)
68	guez	3	END DO
69			DO ij = 1, ip1jm
70	guez	31	pbarvc(ij, l) = pbarvc(ij, l) + pbarv(ij, l)
71	guez	3	END DO
72			END DO
73
74			! selection de la masse instantannee des mailles avant le transport.
75			IF (iadvtr==0) THEN
76	guez	31	CALL scopy(ip1jmp1*llm, masse, 1, massem, 1)
77	guez	3	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	guez	31	CALL groupe(massem, pbaruc, pbarvc, pbarug, pbarvg, wg)
89	guez	3
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	guez	31	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	guez	3	END DO
96	guez	31	CALL scopy(jjm-1, zdp(iip1+iip1), iip1, zdp(iip2), iip1)
97	guez	3	DO ij = iip2, ip1jm
98	guez	31	zdp(ij) = zdp(ij)*dtvr/massem(ij, l)
99	guez	3	END DO
100
101	guez	31	CALL minmax(ip1jm-iip1, zdp(iip2), zdpmin, zdpmax)
102	guez	3
103	guez	31	IF (max(abs(zdpmin), abs(zdpmax))>0.5) THEN
104	guez	3	PRINT *, 'WARNING DP/P l=', l, ' MIN:', zdpmin, ' MAX:', zdpmax
105			END IF
106			END DO
107
108	guez	31	! Advection proprement dite
109	guez	3
110			! Calcul des moyennes basées sur la masse
111	guez	31	CALL massbar(massem, massebx, masseby)
112	guez	3
113			! Appel des sous programmes d'advection
114	guez	31
115	guez	3	DO iq = 1, nqmx
116			IF (iadv(iq)==0) CYCLE
117	guez	31
118	guez	3	! Schema de Van Leer I MUSCL
119	guez	31
120	guez	3	IF (iadv(iq)==10) THEN
121	guez	31	CALL vlsplt(q(1, 1, iq), 2., massem, wg, pbarug, pbarvg, dtvr)
122	guez	3	! Schema "pseudo amont" + test sur humidite specifique
123			! pour la vapeur d'eau. F. Codron
124			ELSE IF (iadv(iq)==14) THEN
125	guez	31	CALL vlspltqs(q(1, 1, 1), 2., massem, wg, pbarug, pbarvg, dtvr, p, &
126			pk, teta)
127	guez	3	! Schema de Frederic Hourdin
128			ELSE IF (iadv(iq)==12) THEN
129			! Pas de temps adaptatif
130	guez	31	CALL adaptdt(iadv(iq), dtbon, n, pbarug, massem)
131	guez	3	IF (n>1) THEN
132	guez	31	WRITE (, ) 'WARNING horizontal dt=', dtbon, 'dtvr=', dtvr, &
133	guez	3	'n=', n
134			END IF
135			DO indice = 1, n
136	guez	31	CALL advn(q(1, 1, iq), massem, wg, pbarug, pbarvg, dtbon, 1)
137	guez	3	END DO
138			ELSE IF (iadv(iq)==13) THEN
139			! Pas de temps adaptatif
140	guez	31	CALL adaptdt(iadv(iq), dtbon, n, pbarug, massem)
141	guez	3	IF (n>1) THEN
142	guez	31	WRITE (, ) 'WARNING horizontal dt=', dtbon, 'dtvr=', dtvr, &
143	guez	3	'n=', n
144			END IF
145			DO indice = 1, n
146	guez	31	CALL advn(q(1, 1, iq), massem, wg, pbarug, pbarvg, dtbon, 2)
147	guez	3	END DO
148			! Schema de pente SLOPES
149			ELSE IF (iadv(iq)==20) THEN
150	guez	31	CALL pentes_ini(q(1, 1, iq), wg, massem, pbarug, pbarvg, 0)
151	guez	3	! Schema de Prather
152			ELSE IF (iadv(iq)==30) THEN
153			! Pas de temps adaptatif
154	guez	31	CALL adaptdt(iadv(iq), dtbon, n, pbarug, massem)
155	guez	3	IF (n>1) THEN
156	guez	31	WRITE (, ) 'WARNING horizontal dt=', dtbon, 'dtvr=', dtvr, &
157	guez	3	'n=', n
158			END IF
159	guez	31	CALL prather(q(1, 1, iq), wg, massem, pbarug, pbarvg, n, dtbon)
160	guez	3	! 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	guez	31	CALL adaptdt(iadv(iq), dtbon, n, pbarug, massem)
165	guez	3	IF (n>1) THEN
166	guez	31	WRITE (, ) 'WARNING horizontal dt=', dtbon, 'dtvr=', dtvr, &
167	guez	3	'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	guez	31	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	guez	3	END DO
178			END DO
179			IF (cflmaxz>=1) THEN
180	guez	31	WRITE (, ) 'WARNING vertical', 'CFLmaxz=', cflmaxz
181	guez	3	END IF
182
183			! Ss-prg interface LMDZ.4->PPM3d
184	guez	31	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	guez	3
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	guez	31	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	guez	3	! Monotonic PPM
196			ELSE IF (iadv(iq)==16) THEN
197			! Ss-prg PPM3d de Lin
198	guez	31	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	guez	3	! Semi Monotonic PPM
202			ELSE IF (iadv(iq)==17) THEN
203			! Ss-prg PPM3d de Lin
204	guez	31	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	guez	3	! Positive Definite PPM
208			ELSE IF (iadv(iq)==18) THEN
209			! Ss-prg PPM3d de Lin
210	guez	31	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	guez	3	END IF
214			END DO
215	guez	31
216	guez	3	! Ss-prg interface PPM3d-LMDZ.4
217	guez	31	CALL interpost(q(1, 1, iq), qppm(1, 1, iq))
218	guez	3	END IF
219			END DO
220
221			! on reinitialise a zero les flux de masse cumules
222			iadvtr = 0
223	guez	23	END IF
224	guez	3
225			END SUBROUTINE advtrac