source: LMDZ6/branches/Amaury_dev/libf/dyn3d/lmdz_advtrac.f90 @ 5218

Last change on this file since 5218 was 5192, checked in by abarral, 3 months ago

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