source: trunk/LMDZ.COMMON/libf/dyn3d/advtrac.F90 @ 1422

Last change on this file since 1422 was 1422, checked in by milmd, 10 years ago

In GENERIC, MARS and COMMON models replace some include files by modules (usefull for decoupling physics with dynamics).

File size: 14.2 KB
RevLine 
[270]1! $Id: advtrac.F90 1549 2011-07-05 08:41:12Z lguez $
[1]2
[270]3SUBROUTINE advtrac(pbaru,pbarv , p,  masse,q,iapptrac,teta, flxw, pk)
4  !     Auteur :  F. Hourdin
5  !
6  !     Modif. P. Le Van     (20/12/97)
7  !            F. Codron     (10/99)
8  !            D. Le Croller (07/2001)
9  !            M.A Filiberti (04/2002)
10  !
[1189]11  USE infotrac, ONLY: nqtot, iadv
12  USE control_mod, ONLY: iapp_tracvl, day_step
[1422]13  USE comconst_mod, ONLY: dtvr
[1]14
15
[270]16  IMPLICIT NONE
17  !
18  include "dimensions.h"
19  include "paramet.h"
20  include "comdissip.h"
21  include "comgeom2.h"
22  include "iniprint.h"
[1]23
[270]24  !-------------------------------------------------------------------
25  !     Arguments
26  !-------------------------------------------------------------------
[1189]27  INTEGER,INTENT(OUT) :: iapptrac
28  REAL,INTENT(IN) :: pbaru(ip1jmp1,llm)
29  REAL,INTENT(IN) :: pbarv(ip1jm,llm)
30  REAL,INTENT(INOUT) :: q(ip1jmp1,llm,nqtot)
31  REAL,INTENT(IN) :: masse(ip1jmp1,llm)
32  REAL,INTENT(IN) :: p( ip1jmp1,llmp1 )
33  REAL,INTENT(IN) :: teta(ip1jmp1,llm)
34  REAL,INTENT(IN) :: pk(ip1jmp1,llm)
35  REAL,INTENT(OUT) :: flxw(ip1jmp1,llm)
36  !-------------------------------------------------------------------
[270]37  !     Ajout PPM
38  !--------------------------------------------------------
39  REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
40  !-------------------------------------------------------------
41  !     Variables locales
42  !-------------------------------------------------------------
[1]43
[270]44  REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
45  REAL massem(ip1jmp1,llm),zdp(ip1jmp1)
46  REAL pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm)
47  REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
48  INTEGER iadvtr
49  INTEGER ij,l,iq,iiq
50  REAL zdpmin, zdpmax
51  EXTERNAL  minmax
52  SAVE iadvtr, massem, pbaruc, pbarvc
53  DATA iadvtr/0/
54  !----------------------------------------------------------
55  !     Rajouts pour PPM
56  !----------------------------------------------------------
57  INTEGER indice,n
58  REAL dtbon ! Pas de temps adaptatif pour que CFL<1
59  REAL CFLmaxz,aaa,bbb ! CFL maximum
60  REAL psppm(iim,jjp1) ! pression  au sol
61  REAL unatppm(iim,jjp1,llm),vnatppm(iim,jjp1,llm)
62  REAL qppm(iim*jjp1,llm,nqtot)
63  REAL fluxwppm(iim,jjp1,llm)
64  REAL apppm(llmp1), bpppm(llmp1)
65  LOGICAL dum,fill
66  DATA fill/.true./
67  DATA dum/.true./
[1]68
[270]69  integer,save :: countcfl=0
70  real cflx(ip1jmp1,llm)
71  real cfly(ip1jm,llm)
72  real cflz(ip1jmp1,llm)
73  real, save :: cflxmax(llm),cflymax(llm),cflzmax(llm)
[1]74
[270]75  IF(iadvtr.EQ.0) THEN
76     CALL initial0(ijp1llm,pbaruc)
77     CALL initial0(ijmllm,pbarvc)
78  ENDIF
[1]79
[270]80  !   accumulation des flux de masse horizontaux
81  DO l=1,llm
82     DO ij = 1,ip1jmp1
83        pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l)
84     ENDDO
85     DO ij = 1,ip1jm
86        pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l)
87     ENDDO
88  ENDDO
[1]89
[270]90  !   selection de la masse instantannee des mailles avant le transport.
91  IF(iadvtr.EQ.0) THEN
[1]92
[270]93     CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
94     !cc         CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
95     !
96  ENDIF
[1]97
[270]98  iadvtr   = iadvtr+1
99  iapptrac = iadvtr
[1]100
101
[270]102  !   Test pour savoir si on advecte a ce pas de temps
103  IF ( iadvtr.EQ.iapp_tracvl ) THEN
[1]104
[270]105     !c   ..  Modif P.Le Van  ( 20/12/97 )  ....
106     !c
[1]107
[270]108     !   traitement des flux de masse avant advection.
109     !     1. calcul de w
110     !     2. groupement des mailles pres du pole.
[1]111
[270]112     CALL groupe( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )
[1]113
[270]114     ! ... Flux de masse diaganostiques traceurs
115     flxw = wg / REAL(iapp_tracvl)
[1]116
[270]117     !  test sur l'eventuelle creation de valeurs negatives de la masse
118     DO l=1,llm-1
119        DO ij = iip2+1,ip1jm
120           zdp(ij) =    pbarug(ij-1,l)   - pbarug(ij,l) &
121                - pbarvg(ij-iip1,l) + pbarvg(ij,l) &
122                +       wg(ij,l+1)  - wg(ij,l)
123        ENDDO
124        CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
125        DO ij = iip2,ip1jm
126           zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
127        ENDDO
[1]128
129
[270]130        CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
[1]131
[270]132        IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN
133           PRINT*,'WARNING DP/P l=',l,'  MIN:',zdpmin, &
134                '   MAX:', zdpmax
135        ENDIF
[1]136
[270]137     ENDDO
[1]138
139
[270]140     !-------------------------------------------------------------------
141     ! Calcul des criteres CFL en X, Y et Z
142     !-------------------------------------------------------------------
[1]143
[270]144     if (countcfl == 0. ) then
145        cflxmax(:)=0.
146        cflymax(:)=0.
147        cflzmax(:)=0.
148     endif
[1]149
[270]150     countcfl=countcfl+iapp_tracvl
151     cflx(:,:)=0.
152     cfly(:,:)=0.
153     cflz(:,:)=0.
154     do l=1,llm
155        do ij=iip2,ip1jm-1
156           if (pbarug(ij,l)>=0.) then
157              cflx(ij,l)=pbarug(ij,l)*dtvr/masse(ij,l)
158           else
159              cflx(ij,l)=-pbarug(ij,l)*dtvr/masse(ij+1,l)
160           endif
161        enddo
162     enddo
163     do l=1,llm
164        do ij=iip2,ip1jm-1,iip1
165           cflx(ij+iip1,l)=cflx(ij,l)
166        enddo
167     enddo
[1]168
[270]169     do l=1,llm
170        do ij=1,ip1jm
171           if (pbarvg(ij,l)>=0.) then
172              cfly(ij,l)=pbarvg(ij,l)*dtvr/masse(ij,l)
173           else
174              cfly(ij,l)=-pbarvg(ij,l)*dtvr/masse(ij+iip1,l)
175           endif
176        enddo
177     enddo
[1]178
[270]179     do l=2,llm
180        do ij=1,ip1jm
181           if (wg(ij,l)>=0.) then
182              cflz(ij,l)=wg(ij,l)*dtvr/masse(ij,l)
183           else
184              cflz(ij,l)=-wg(ij,l)*dtvr/masse(ij,l-1)
185           endif
186        enddo
187     enddo
188
189     do l=1,llm
190        cflxmax(l)=max(cflxmax(l),maxval(cflx(:,l)))
191        cflymax(l)=max(cflymax(l),maxval(cfly(:,l)))
192        cflzmax(l)=max(cflzmax(l),maxval(cflz(:,l)))
193     enddo
194
[1]195!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[270]196     ! Par defaut, on sort le diagnostic des CFL tous les jours.
197     ! Si on veut le sortir a chaque pas d'advection en cas de plantage
198     !     if (countcfl==iapp_tracvl) then
[1]199!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[270]200     if (countcfl==day_step) then
201        do l=1,llm
[1019]202           write(lunout,*) 'L, CFL[xyz]max:', l, cflxmax(l), cflymax(l), &
203                cflzmax(l)
[270]204        enddo
205        countcfl=0
206     endif
[1]207
[270]208     !-------------------------------------------------------------------
209     !   Advection proprement dite (Modification Le Croller (07/2001)
210     !-------------------------------------------------------------------
[1]211
[270]212     !----------------------------------------------------
213     !        Calcul des moyennes basées sur la masse
214     !----------------------------------------------------
215     call massbar(massem,massebx,masseby)         
216
217     !-----------------------------------------------------------
218     !     Appel des sous programmes d'advection
219     !-----------------------------------------------------------
220     do iq=1,nqtot
221        !        call clock(t_initial)
[1]222        if(iadv(iq) == 0) cycle
[270]223        !   ----------------------------------------------------------------
224        !   Schema de Van Leer I MUSCL
225        !   ----------------------------------------------------------------
[1]226        if(iadv(iq).eq.10) THEN
[270]227           call vlsplt(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
[1]228
[270]229           !   ----------------------------------------------------------------
230           !   Schema "pseudo amont" + test sur humidite specifique
231           !    pour la vapeur d'eau. F. Codron
232           !   ----------------------------------------------------------------
[1]233        else if(iadv(iq).eq.14) then
[270]234           !
235           CALL vlspltqs( q(1,1,1), 2., massem, wg , &
236                pbarug,pbarvg,dtvr,p,pk,teta )
237           !   ----------------------------------------------------------------
238           !   Schema de Frederic Hourdin
239           !   ----------------------------------------------------------------
[1]240        else if(iadv(iq).eq.12) then
[270]241           !            Pas de temps adaptatif
[1]242           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
243           if (n.GT.1) then
[270]244              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
245                   dtvr,'n=',n
[1]246           endif
247           do indice=1,n
[270]248              call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
[1]249           end do
250        else if(iadv(iq).eq.13) then
[270]251           !            Pas de temps adaptatif
[1]252           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
253           if (n.GT.1) then
[270]254              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
255                   dtvr,'n=',n
[1]256           endif
[270]257           do indice=1,n
258              call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
259           end do
260           !   ----------------------------------------------------------------
261           !   Schema de pente SLOPES
262           !   ----------------------------------------------------------------
[1]263        else if (iadv(iq).eq.20) then
[270]264           call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
[1]265
[270]266           !   ----------------------------------------------------------------
267           !   Schema de Prather
268           !   ----------------------------------------------------------------
[1]269        else if (iadv(iq).eq.30) then
[270]270           !            Pas de temps adaptatif
[1]271           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
272           if (n.GT.1) then
[270]273              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
274                   dtvr,'n=',n
[1]275           endif
[270]276           call  prather(q(1,1,iq),wg,massem,pbarug,pbarvg, &
277                n,dtbon)
[1]278
[270]279           !   ----------------------------------------------------------------
280           !   Schemas PPM Lin et Rood
281           !   ----------------------------------------------------------------
282        else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND. &
283             iadv(iq).LE.18)) then
[1]284
[270]285           !        Test sur le flux horizontal
286           !        Pas de temps adaptatif
287           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
288           if (n.GT.1) then
289              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
290                   dtvr,'n=',n
291           endif
292           !        Test sur le flux vertical
293           CFLmaxz=0.
294           do l=2,llm
295              do ij=iip2,ip1jm
296                 aaa=wg(ij,l)*dtvr/massem(ij,l)
297                 CFLmaxz=max(CFLmaxz,aaa)
298                 bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
299                 CFLmaxz=max(CFLmaxz,bbb)
300              enddo
[1]301           enddo
[270]302           if (CFLmaxz.GE.1) then
303              write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
304           endif
[1]305
[270]306           !-----------------------------------------------------------
307           !        Ss-prg interface LMDZ.4->PPM3d
308           !-----------------------------------------------------------
[1]309
[270]310           call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, &
311                apppm,bpppm,massebx,masseby,pbarug,pbarvg, &
312                unatppm,vnatppm,psppm)
[1]313
[270]314           do indice=1,n
315              !----------------------------------------------------------------
316              !                         VL (version PPM) horiz. et PPM vert.
317              !----------------------------------------------------------------
318              if (iadv(iq).eq.11) then
319                 !                  Ss-prg PPM3d de Lin
320                 call ppm3d(1,qppm(1,1,iq), &
321                      psppm,psppm, &
322                      unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1, &
323                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
324                      fill,dum,220.)
[1]325
[270]326                 !-------------------------------------------------------------
327                 !                           Monotonic PPM
328                 !-------------------------------------------------------------
329              else if (iadv(iq).eq.16) then
330                 !                  Ss-prg PPM3d de Lin
331                 call ppm3d(1,qppm(1,1,iq), &
332                      psppm,psppm, &
333                      unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1, &
334                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
335                      fill,dum,220.)
336                 !-------------------------------------------------------------
[1]337
[270]338                 !-------------------------------------------------------------
339                 !                           Semi Monotonic PPM
340                 !-------------------------------------------------------------
341              else if (iadv(iq).eq.17) then
342                 !                  Ss-prg PPM3d de Lin
343                 call ppm3d(1,qppm(1,1,iq), &
344                      psppm,psppm, &
345                      unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1, &
346                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
347                      fill,dum,220.)
348                 !-------------------------------------------------------------
[1]349
[270]350                 !-------------------------------------------------------------
351                 !                         Positive Definite PPM
352                 !-------------------------------------------------------------
353              else if (iadv(iq).eq.18) then
354                 !                  Ss-prg PPM3d de Lin
355                 call ppm3d(1,qppm(1,1,iq), &
356                      psppm,psppm, &
357                      unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1, &
358                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
359                      fill,dum,220.)
360                 !-------------------------------------------------------------
361              endif
362           enddo
363           !-----------------------------------------------------------------
364           !               Ss-prg interface PPM3d-LMDZ.4
365           !-----------------------------------------------------------------
366           call interpost(q(1,1,iq),qppm(1,1,iq))
367        endif
368        !----------------------------------------------------------------------
[1]369
[270]370        !-----------------------------------------------------------------
371        ! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
372        ! et Nord j=1
373        !-----------------------------------------------------------------
[1]374
[270]375        !                  call traceurpole(q(1,1,iq),massem)
[1]376
[270]377        ! calcul du temps cpu pour un schema donne
[1]378
[270]379        !                  call clock(t_final)
380        !ym                  tps_cpu=t_final-t_initial
381        !ym                  cpuadv(iq)=cpuadv(iq)+tps_cpu
[1]382
[270]383     end DO
[1]384
385
[270]386     !------------------------------------------------------------------
387     !   on reinitialise a zero les flux de masse cumules
388     !---------------------------------------------------
389     iadvtr=0
[1]390
[270]391  ENDIF ! if iadvtr.EQ.iapp_tracvl
[1]392
[270]393END SUBROUTINE advtrac
Note: See TracBrowser for help on using the repository browser.