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

Diff of /trunk/dyn3d/advtrac.f

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

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