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 dimensions, ONLY : jjm, llm, nqmx
    use groupe_m, only: groupe
    USE iniadvtrac_m, ONLY : iadv
    use massbar_m, only: massbar
    USE paramet_m, ONLY : iip1, iip2, ijmllm, ijp1llm, ip1jm, ip1jmp1, llmp1
    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

    INTEGER indice, n
    ! Pas de temps adaptatif pour que CFL < 1 
    REAL dtbon

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

    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 (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 (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)
          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	265	USE dimensions, ONLY : jjm, llm, nqmx
15	guez	207	use groupe_m, only: groupe
16	guez	71	USE iniadvtrac_m, ONLY : iadv
17	guez	91	use massbar_m, only: massbar
18	guez	178	USE paramet_m, ONLY : iip1, iip2, ijmllm, ijp1llm, ip1jm, ip1jmp1, llmp1
19	guez	157	use vlsplt_m, only: vlsplt
20	guez	108	use vlspltqs_m, only: vlspltqs
21	guez	31
22	guez	71	REAL, intent(in):: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
23			REAL, intent(in):: p(ip1jmp1, llmp1)
24			real, intent(in):: masse(ip1jmp1, llm)
25	guez	40	REAL, intent(inout):: q(ip1jmp1, llm, nqmx)
26	guez	71	INTEGER, intent(out):: iapptrac
27	guez	44	real, intent(in):: teta(ip1jmp1, llm)
28	guez	71	REAL, intent(in):: pk(ip1jmp1, llm)
29	guez	3
30	guez	40	! Variables locales
31	guez	3
32	guez	40	REAL massebx(ip1jmp1, llm), masseby(ip1jm, llm)
33			REAL, save:: pbaruc(ip1jmp1, llm), pbarvc(ip1jm, llm)
34			REAL, save:: massem(ip1jmp1, llm)
35			real zdp(ip1jmp1)
36			REAL pbarug(ip1jmp1, llm), pbarvg(ip1jm, llm), wg(ip1jmp1, llm)
37	guez	31
38	guez	40	INTEGER:: iadvtr = 0
39			INTEGER ij, l, iq
40			REAL zdpmin, zdpmax
41			EXTERNAL minmax
42	guez	31
43	guez	40	INTEGER indice, n
44			! Pas de temps adaptatif pour que CFL < 1
45			REAL dtbon
46	guez	3
47	guez	40	!-----------------------------------------------------------
48	guez	3
49	guez	40	IF (iadvtr==0) THEN
50			CALL initial0(ijp1llm, pbaruc)
51			CALL initial0(ijmllm, pbarvc)
52			END IF
53	guez	3
54	guez	40	! accumulation des flux de masse horizontaux
55			DO l = 1, llm
56			DO ij = 1, ip1jmp1
57			pbaruc(ij, l) = pbaruc(ij, l) + pbaru(ij, l)
58			END DO
59			DO ij = 1, ip1jm
60			pbarvc(ij, l) = pbarvc(ij, l) + pbarv(ij, l)
61			END DO
62			END DO
63	guez	3
64	guez	40	! selection de la masse instantannee des mailles avant le transport.
65	guez	71	IF (iadvtr==0) massem = masse
66	guez	3
67	guez	40	iadvtr = iadvtr + 1
68			iapptrac = iadvtr
69	guez	3
70	guez	40	! Test pour savoir si on advecte a ce pas de temps
71	guez	71	IF (iadvtr == iapp_tracvl) THEN
72	guez	40	! traitement des flux de masse avant advection.
73			! 1. calcul de w
74			! 2. groupement des mailles pres du pole.
75	guez	3
76	guez	150	CALL groupe(pbaruc, pbarvc, pbarug, pbarvg, wg)
77	guez	3
78	guez	40	! test sur l'eventuelle creation de valeurs negatives de la masse
79			DO l = 1, llm - 1
80			DO ij = iip2 + 1, ip1jm
81			zdp(ij) = pbarug(ij-1, l) - pbarug(ij, l) - pbarvg(ij-iip1, l) + &
82			pbarvg(ij, l) + wg(ij, l+1) - wg(ij, l)
83			END DO
84			CALL scopy(jjm-1, zdp(iip1+iip1), iip1, zdp(iip2), iip1)
85			DO ij = iip2, ip1jm
86			zdp(ij) = zdp(ij)*dtvr/massem(ij, l)
87			END DO
88	guez	3
89	guez	40	CALL minmax(ip1jm-iip1, zdp(iip2), zdpmin, zdpmax)
90	guez	3
91	guez	40	IF (max(abs(zdpmin), abs(zdpmax))>0.5) THEN
92			PRINT *, 'WARNING DP/P l=', l, ' MIN:', zdpmin, ' MAX:', zdpmax
93			END IF
94			END DO
95	guez	3
96	guez	40	! Advection proprement dite
97	guez	3
98	guez	150	! Calcul des moyennes bas\'ees sur la masse
99	guez	40	CALL massbar(massem, massebx, masseby)
100	guez	31
101	guez	40	! Appel des sous programmes d'advection
102	guez	31
103	guez	40	DO iq = 1, nqmx
104	guez	157	select case (iadv(iq))
105			case (10)
106			! Schema de Van Leer I MUSCL
107	guez	40	CALL vlsplt(q(:, :, iq), 2., massem, wg, pbarug, pbarvg, dtvr)
108	guez	177	case (12)
109	guez	40	! Schema de Frederic Hourdin
110			! Pas de temps adaptatif
111	guez	104	CALL adaptdt(dtbon, n, pbarug, massem)
112	guez	40	IF (n>1) THEN
113			WRITE (, ) 'WARNING horizontal dt=', dtbon, 'dtvr=', dtvr, &
114			'n=', n
115			END IF
116			DO indice = 1, n
117			CALL advn(q(1, 1, iq), massem, wg, pbarug, pbarvg, dtbon, 1)
118			END DO
119	guez	157	case (13)
120	guez	40	! Pas de temps adaptatif
121	guez	104	CALL adaptdt(dtbon, n, pbarug, massem)
122	guez	40	IF (n>1) THEN
123			WRITE (, ) 'WARNING horizontal dt=', dtbon, 'dtvr=', dtvr, &
124			'n=', n
125			END IF
126			DO indice = 1, n
127			CALL advn(q(1, 1, iq), massem, wg, pbarug, pbarvg, dtbon, 2)
128			END DO
129	guez	177	case (14)
130			! Schema "pseudo amont" + test sur humidite specifique
131			! pour la vapeur d'eau. F. Codron
132			CALL vlspltqs(q(1, 1, 1), 2., massem, wg, pbarug, pbarvg, dtvr, &
133			p, pk, teta)
134	guez	157	END select
135	guez	40	END DO
136	guez	3
137	guez	40	! on reinitialise a zero les flux de masse cumules
138			iadvtr = 0
139			END IF
140
141			END SUBROUTINE advtrac
142
143			end module advtrac_m