/[lmdze]/trunk/dyn3d/advtrac.f
ViewVC logotype

Diff of /trunk/dyn3d/advtrac.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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