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

Diff of /trunk/Sources/dyn3d/advtrac.f

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

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