source: LMDZ6/trunk/libf/dyn3d/advtrac.f90 @ 5273

Last change on this file since 5273 was 5272, checked in by abarral, 2 days ago

Turn paramet.h into a module

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.4 KB
Line 
1! $Id: advtrac.f90 5272 2024-10-24 15:53:15Z abarral $
2
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   !
11   USE infotrac,     ONLY: nqtot, tracers, isoCheck
12   USE control_mod,  ONLY: iapp_tracvl, day_step
13   USE comconst_mod, ONLY: dtvr
14   USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_DEBUGIO
15   USE strings_mod, ONLY: int2str
16   USE dimensions_mod, ONLY: iim, jjm, llm, ndm
17USE paramet_mod_h, ONLY: iip1, iip2, iip3, jjp1, llmp1, llmp2, llmm1, kftd, ip1jm, ip1jmp1, &
18          ip1jmi1, ijp1llm, ijmllm, mvar, jcfil, jcfllm
19IMPLICIT NONE
20   !
21
22
23   include "comdissip.h"
24   include "comgeom2.h"
25   include "description.h"
26   include "iniprint.h"
27
28   !---------------------------------------------------------------------------
29   !     Arguments
30   !---------------------------------------------------------------------------
31   INTEGER, INTENT(OUT) :: iapptrac
32   REAL, INTENT(IN) :: pbaru(ip1jmp1,llm)
33   REAL, INTENT(IN) :: pbarv(ip1jm,  llm)
34   REAL, INTENT(INOUT) ::  q(ip1jmp1,llm,nqtot)
35   REAL, INTENT(IN) :: masse(ip1jmp1,llm)
36   REAL, INTENT(IN) ::     p(ip1jmp1,llmp1 )
37   REAL, INTENT(IN) ::  teta(ip1jmp1,llm)
38   REAL, INTENT(IN) ::    pk(ip1jmp1,llm)
39   REAL, INTENT(OUT) :: flxw(ip1jmp1,llm)
40   !---------------------------------------------------------------------------
41   !     Ajout PPM
42   !---------------------------------------------------------------------------
43   REAL :: massebx(ip1jmp1,llm), masseby(ip1jm,llm)
44   !---------------------------------------------------------------------------
45   !     Variables locales
46   !---------------------------------------------------------------------------
47   INTEGER :: ij, l, iq, iadv
48!   REAL(KIND=KIND(1.d0)) :: t_initial, t_final, tps_cpu
49   REAL :: zdp(ip1jmp1), zdpmin, zdpmax
50   INTEGER, SAVE :: iadvtr=0
51   REAL, DIMENSION(ip1jmp1,llm) :: pbaruc, pbarug, massem, wg
52   REAL, DIMENSION(ip1jm,  llm) :: pbarvc, pbarvg
53   EXTERNAL  minmax
54   SAVE massem, pbaruc, pbarvc
55   !---------------------------------------------------------------------------
56   !     Rajouts pour PPM
57   !---------------------------------------------------------------------------
58   INTEGER indice, n
59   REAL :: dtbon                       ! Pas de temps adaptatif pour que CFL<1
60   REAL :: CFLmaxz, aaa, bbb           ! CFL maximum
61   REAL, DIMENSION(iim,jjp1,llm) :: unatppm, vnatppm, fluxwppm
62   REAL ::    qppm(iim*jjp1,llm,nqtot)
63   REAL ::   psppm(iim,jjp1)           ! pression  au sol
64   REAL, DIMENSION(llmp1) :: apppm, bpppm
65   LOGICAL, SAVE :: dum=.TRUE., fill=.TRUE.
66
67   INTEGER, SAVE :: countcfl=0
68   REAL, DIMENSION(ip1jmp1,llm) :: cflx, cflz
69   REAL, DIMENSION(ip1jm  ,llm) :: cfly
70   REAL, DIMENSION(llm), SAVE :: cflxmax, cflymax, cflzmax
71
72   IF(iadvtr == 0) THEN
73      pbaruc(:,:)=0
74      pbarvc(:,:)=0
75   END IF
76
77   !--- Accumulation des flux de masse horizontaux
78   DO l=1,llm
79      DO ij = 1,ip1jmp1
80         pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l)
81      END DO
82      DO ij = 1,ip1jm
83         pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l)
84      END DO
85   END DO
86
87   !--- Selection de la masse instantannee des mailles avant le transport.
88   IF(iadvtr == 0) THEN
89     CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
90   ! CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
91   END IF
92
93   iadvtr   = iadvtr+1
94   iapptrac = iadvtr
95
96   !--- Test pour savoir si on advecte a ce pas de temps
97   IF(iadvtr /= iapp_tracvl) RETURN
98
99   !   ..  Modif P.Le Van  ( 20/12/97 )  ....
100   !
101   !   traitement des flux de masse avant advection.
102   !       1. calcul de w
103   !       2. groupement des mailles pres du pole.
104
105   CALL groupe(massem, pbaruc, pbarvc, pbarug, pbarvg, wg)
106
107   !--- Flux de masse diaganostiques traceurs
108   flxw = wg / REAL(iapp_tracvl)
109
110   !--- Test sur l'eventuelle creation de valeurs negatives de la masse
111   DO l=1,llm-1
112      DO ij = iip2+1,ip1jm
113         zdp(ij) = pbarug(ij-1,l)    - pbarug(ij,l) &
114                 - pbarvg(ij-iip1,l) + pbarvg(ij,l) &
115                 +     wg(ij,l+1)    -     wg(ij,l)
116      END DO
117! ym  ---> pourquoi jjm-1 et non jjm ? a cause du pole ?
118      CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
119      DO ij = iip2,ip1jm
120         zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
121      END DO
122
123      CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
124
125      IF(MAX(ABS(zdpmin),ABS(zdpmax)) > 0.5) &
126         WRITE(*,*)'WARNING DP/P l=',l,'  MIN:',zdpmin,' MAX:', zdpmax
127
128   END DO
129
130   !-------------------------------------------------------------------------
131   ! Calcul des criteres CFL en X, Y et Z
132   !-------------------------------------------------------------------------
133   IF(countcfl == 0. ) then
134      cflxmax(:)=0.
135      cflymax(:)=0.
136      cflzmax(:)=0.
137   END IF
138
139   countcfl=countcfl+iapp_tracvl
140   cflx(:,:)=0.
141   cfly(:,:)=0.
142   cflz(:,:)=0.
143   DO l=1,llm
144      DO ij=iip2,ip1jm-1
145         IF(pbarug(ij,l)>=0.) then
146            cflx(ij,l)=pbarug(ij,l)*dtvr/masse(ij,l)
147         ELSE
148            cflx(ij,l)=-pbarug(ij,l)*dtvr/masse(ij+1,l)
149         END IF
150      END DO
151   END DO
152
153   DO l=1,llm
154      DO ij=iip2,ip1jm-1,iip1
155         cflx(ij+iip1,l)=cflx(ij,l)
156      END DO
157   END DO
158
159   DO l=1,llm
160      DO ij=1,ip1jm
161         IF(pbarvg(ij,l)>=0.) then
162            cfly(ij,l)=pbarvg(ij,l)*dtvr/masse(ij,l)
163         ELSE
164            cfly(ij,l)=-pbarvg(ij,l)*dtvr/masse(ij+iip1,l)
165         END IF
166      END DO
167   END DO
168
169   DO l=2,llm
170      DO ij=1,ip1jm
171         IF(wg(ij,l) >= 0.) THEN
172            cflz(ij,l)=wg(ij,l)*dtvr/masse(ij,l)
173         ELSE
174            cflz(ij,l)=-wg(ij,l)*dtvr/masse(ij,l-1)
175         END IF
176      END DO
177   END DO
178
179   DO l=1,llm
180      cflxmax(l)=max(cflxmax(l),maxval(cflx(:,l)))
181      cflymax(l)=max(cflymax(l),maxval(cfly(:,l)))
182      cflzmax(l)=max(cflzmax(l),maxval(cflz(:,l)))
183   END DO
184
185   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
186   ! Par defaut, on sort le diagnostic des CFL tous les jours.
187   ! Si on veut le sortir a chaque pas d'advection en cas de plantage
188   !       IF(countcfl==iapp_tracvl) then
189   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
190   IF(countcfl==day_step) then
191      DO l=1,llm
192         WRITE(lunout,*) 'L, CFL[xyz]max:', l, cflxmax(l), cflymax(l), cflzmax(l)
193      END DO
194      countcfl=0
195   END IF
196
197   !---------------------------------------------------------------------------
198   !   Advection proprement dite (Modification Le Croller (07/2001)
199   !---------------------------------------------------------------------------
200
201   !---------------------------------------------------------------------------
202   !   Calcul des moyennes basees sur la masse
203   !---------------------------------------------------------------------------
204   CALL massbar(massem,massebx,masseby)
205
206IF (CPPKEY_DEBUGIO) THEN
207   CALL WriteField_u('massem',massem)
208   CALL WriteField_u('wg',wg)
209   CALL WriteField_u('pbarug',pbarug)
210   CALL WriteField_v('pbarvg',pbarvg)
211   CALL WriteField_u('p_tmp',p)
212   CALL WriteField_u('pk_tmp',pk)
213   CALL WriteField_u('teta_tmp',teta)
214   DO iq=1,nqtot
215      CALL WriteField_u('q_adv'//trim(int2str(iq)),q(:,:,iq))
216   END DO
217END IF
218
219   IF(isoCheck) WRITE(*,*) 'advtrac 227'
220   CALL check_isotopes_seq(q,ip1jmp1,'advtrac 162')
221
222   !-------------------------------------------------------------------------
223   !       Appel des sous programmes d'advection
224   !-------------------------------------------------------------------------
225   DO iq = 1, nqtot
226!     CALL clock(t_initial)
227      IF(tracers(iq)%parent /= 'air') CYCLE
228      iadv = tracers(iq)%iadv
229      !-----------------------------------------------------------------------
230      SELECT CASE(iadv)
231      !-----------------------------------------------------------------------
232         CASE(0); CYCLE
233         !--------------------------------------------------------------------
234         CASE(10)  !--- Schema de Van Leer I MUSCL
235         !--------------------------------------------------------------------
236!           WRITE(*,*) 'advtrac 239: iq,q(1721,19,:)=',iq,q(1721,19,:)     
237            CALL vlsplt(q,2.,massem,wg,pbarug,pbarvg,dtvr,iq)
238
239         !--------------------------------------------------------------------
240         CASE(14)  !--- Schema "pseuDO amont" + test sur humidite specifique
241                   !--- pour la vapeur d'eau. F. Codron
242         !--------------------------------------------------------------------
243!           WRITE(*,*) 'advtrac 248: iq,q(1721,19,:)=',iq,q(1721,19,:)
244            CALL vlspltqs(q,2.,massem,wg,pbarug,pbarvg,dtvr,p,pk,teta,iq)
245
246         !--------------------------------------------------------------------
247         CASE(12)  !--- Schema de Frederic Hourdin
248         !--------------------------------------------------------------------
249            CALL adaptdt(iadv,dtbon,n,pbarug,massem)   ! pas de temps adaptatif
250            IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n
251            DO indice=1,n
252              CALL advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
253            END DO
254
255         !--------------------------------------------------------------------
256         CASE(13)  !--- Pas de temps adaptatif
257         !--------------------------------------------------------------------
258            CALL adaptdt(iadv,dtbon,n,pbarug,massem)
259            IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n
260            DO indice=1,n
261               CALL advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
262            END DO
263
264         !--------------------------------------------------------------------
265         CASE(20)  !--- Schema de pente SLOPES
266         !--------------------------------------------------------------------
267            CALL pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
268
269         !--------------------------------------------------------------------
270         CASE(30)  !--- Schema de Prather
271         !--------------------------------------------------------------------
272            ! Pas de temps adaptatif
273            CALL adaptdt(iadv,dtbon,n,pbarug,massem)
274            IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n
275            CALL prather(q(1,1,iq),wg,massem,pbarug,pbarvg,n,dtbon)
276
277         !--------------------------------------------------------------------
278         CASE(11,16,17,18)   !--- Schemas PPM Lin et Rood
279         !--------------------------------------------------------------------
280            ! Test sur le flux horizontal
281            CALL adaptdt(iadv,dtbon,n,pbarug,massem)   ! pas de temps adaptatif
282            IF(n > 1) WRITE(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',dtvr,'n=',n
283            ! Test sur le flux vertical
284            CFLmaxz=0.
285            DO l=2,llm
286               DO ij=iip2,ip1jm
287                  aaa=wg(ij,l)*dtvr/massem(ij,l)
288                  CFLmaxz=max(CFLmaxz,aaa)
289                  bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
290                  CFLmaxz=max(CFLmaxz,bbb)
291               END DO
292            END DO
293            IF(CFLmaxz.GE.1) WRITE(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
294            !----------------------------------------------------------------
295            !     Ss-prg interface LMDZ.4->PPM3d (ss-prg de Lin)
296            !----------------------------------------------------------------
297            CALL interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, &
298                 apppm,bpppm,massebx,masseby,pbarug,pbarvg, &
299                 unatppm,vnatppm,psppm)
300
301            !----------------------------------------------------------------
302            DO indice=1,n     !--- VL (version PPM) horiz. et PPM vert.
303            !----------------------------------------------------------------
304               SELECT CASE(iadv)
305                  !----------------------------------------------------------
306                  CASE(11)
307                  !----------------------------------------------------------
308                     CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, &
309                                2,2,2,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,fill,dum,220.)
310                  !----------------------------------------------------------
311                  CASE(16) !--- Monotonic PPM
312                  !----------------------------------------------------------
313                     CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, &
314                                3,3,3,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,fill,dum,220.)
315                  !----------------------------------------------------------
316                  CASE(17) !--- Semi monotonic PPM
317                  !----------------------------------------------------------
318                     CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, &
319                                4,4,4,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, fill,dum,220.)
320                  !----------------------------------------------------------
321                  CASE(18) !--- Positive Definite PPM
322                  !----------------------------------------------------------
323                     CALL ppm3d(1,qppm(1,1,iq),psppm,psppm,unatppm,vnatppm,fluxwppm,dtbon, &
324                                5,5,5,1,iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,fill,dum,220.)
325               END SELECT
326            !----------------------------------------------------------------
327            END DO
328            !----------------------------------------------------------------
329            !     Ss-prg interface PPM3d-LMDZ.4
330            !----------------------------------------------------------------
331            CALL interpost(q(1,1,iq),qppm(1,1,iq))
332      !----------------------------------------------------------------------
333      END SELECT
334      !----------------------------------------------------------------------
335
336      !----------------------------------------------------------------------
337      ! On impose une seule valeur du traceur au pole Sud j=jjm+1=jjp1 et Nord j=1
338      !----------------------------------------------------------------------
339      !  CALL traceurpole(q(1,1,iq),massem)
340
341      !--- Calcul du temps cpu pour un schema donne
342      !  CALL clock(t_final)
343      !ym  tps_cpu=t_final-t_initial
344      !ym  cpuadv(iq)=cpuadv(iq)+tps_cpu
345
346   END DO
347
348   IF(isoCheck) WRITE(*,*) 'advtrac 402'
349   CALL check_isotopes_seq(q,ip1jmp1,'advtrac 397')
350
351   !-------------------------------------------------------------------------
352   !   on reinitialise a zero les flux de masse cumules
353   !-------------------------------------------------------------------------
354   iadvtr=0
355
356END SUBROUTINE advtrac
Note: See TracBrowser for help on using the repository browser.