source: LMDZ6/branches/Amaury_dev/libf/dyn3dmem/advtrac_loc.f90 @ 5225

Last change on this file since 5225 was 5195, checked in by abarral, 2 months ago

Correct r5192, some lmdz_description cases were missing
  • 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
File size: 12.4 KB
RevLine 
[4056]1SUBROUTINE advtrac_loc(pbarug, pbarvg, wg, p, massem, q, teta, pk)
[5101]2  !     Auteur :  F. Hourdin
[5099]3
[5101]4  !     Modif. P. Le Van     (20/12/97)
5  !            F. Codron     (10/99)
6  !            D. Le Croller (07/2001)
7  !            M.A Filiberti (04/2002)
[5099]8
[5182]9  USE lmdz_infotrac, ONLY: nqtot, tracers
[5101]10  USE control_mod, ONLY: iapp_tracvl, day_step, planet_type
11  USE comconst_mod, ONLY: dtvr
12  USE parallel_lmdz
13  USE Write_Field_loc
14  USE Bands
15  USE mod_hallo
[5117]16  USE lmdz_vampir
[5101]17  USE times
18  USE advtrac_mod, ONLY: finmasse
19  USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_DEBUGIO
[5117]20  USE lmdz_strings, ONLY: int2str
[5116]21  USE lmdz_libmath, ONLY: minmax
[5134]22  USE lmdz_comdissip, ONLY: coefdis, tetavel, tetatemp, gamdissip, niterdis
[5136]23  USE lmdz_comgeom2
[1632]24
[5159]25USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
26  USE lmdz_paramet
[5101]27  IMPLICIT NONE
[5099]28
[1632]29
[5159]30
31
[5101]32  !---------------------------------------------------------------------------
33  !     Arguments
34  !---------------------------------------------------------------------------
35  REAL, INTENT(IN) :: pbarug(ijb_u:ije_u, llm)
36  REAL, INTENT(IN) :: pbarvg(ijb_v:ije_v, llm)
37  REAL, INTENT(IN) :: wg(ijb_u:ije_u, llm)
38  REAL, INTENT(IN) :: p(ijb_u:ije_u, llmp1)
39  REAL, INTENT(IN) :: massem(ijb_u:ije_u, llm)
40  REAL, INTENT(INOUT) :: q(ijb_u:ije_u, llm, nqtot)
41  REAL, INTENT(IN) :: teta(ijb_u:ije_u, llm)
42  REAL, INTENT(IN) :: pk(ijb_u:ije_u, llm)
43  !---------------------------------------------------------------------------
44  !     Ajout PPM
45  !---------------------------------------------------------------------------
46  REAL :: massebx(ijb_u:ije_u, llm), masseby(ijb_v:ije_v, llm)
47  !---------------------------------------------------------------------------
48  !     Variables locales
49  !---------------------------------------------------------------------------
50  INTEGER :: ij, l, iq, iadv
51  REAL(KIND = KIND(1.d0)) :: t_initial, t_final, tps_cpu
52  REAL :: zdp(ijb_u:ije_u), zdpmin, zdpmax
53  INTEGER, SAVE :: iadvtr = 0
54  !$OMP THREADPRIVATE(iadvtr)
[1632]55
[5101]56  !---------------------------------------------------------------------------
57  !     Rajouts pour PPM
58  !---------------------------------------------------------------------------
59  INTEGER :: indice, n
60  REAL :: dtbon                       ! Pas de temps adaptatif pour que CFL<1
61  REAL :: CFLmaxz, aaa, bbb           ! CFL maximum
62  REAL, DIMENSION(iim, jjb_u:jje_u, llm) :: unatppm, vnatppm, fluxwppm
63  REAL :: qppm(iim * jjnb_u, llm, nqtot)
64  REAL :: psppm(iim, jjb_u:jje_u)    ! pression  au sol
65  REAL, DIMENSION(llmp1) :: apppm, bpppm
66  LOGICAL, SAVE :: dum = .TRUE., fill = .TRUE.
67  INTEGER :: ijb, ije, ijbu, ijbv, ijeu, ijev, j
68  TYPE(Request), SAVE :: testRequest
69  !$OMP THREADPRIVATE(testRequest)
[1632]70
[5101]71  ! Test sur l'eventuelle creation de valeurs negatives de la masse
72  ijb = ij_begin; IF(pole_nord) ijb = ij_begin + iip1
73  ije = ij_end;   IF(pole_sud)  ije = ij_end - iip1
[1632]74
[5101]75  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
76  DO l = 1, llm - 1
77    DO ij = ijb + 1, ije
78      zdp(ij) = pbarug(ij - 1, l) - pbarug(ij, l) &
79              - pbarvg(ij - iip1, l) + pbarvg(ij, l) &
80              + wg(ij, l + 1) - wg(ij, l)
81    END DO
82    ! ym  ---> pourquoi jjm-1 et non jjm ? a cause du pole ?
83    !     CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
84    DO ij = ijb, ije - iip1 + 1, iip1
85      zdp(ij) = zdp(ij + iip1 - 1)
86    END DO
87    DO ij = ijb, ije
88      zdp(ij) = zdp(ij) * dtvr / massem(ij, l)
89    END DO
90    !     CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
91    ! ym ---> eventuellement a revoir
92    CALL minmax(ije - ijb + 1, zdp(ijb), zdpmin, zdpmax)
93    IF(MAX(ABS(zdpmin), ABS(zdpmax)) >0.5) &
94            WRITE(*, *)'WARNING DP/P l=', l, '  MIN:', zdpmin, '   MAX:', zdpmax
95  END DO
96  !$OMP END DO NOWAIT
[1632]97
[5101]98  !---------------------------------------------------------------------------
99  !   Advection proprement dite (Modification Le Croller (07/2001)
100  !---------------------------------------------------------------------------
[1632]101
[5101]102  !---------------------------------------------------------------------------
103  !   Calcul des moyennes basees sur la masse
104  !---------------------------------------------------------------------------
105  !ym   CALL massbar_p(massem,massebx,masseby)
106  !ym   ----> Normalement, inutile pour les schemas classiques
107  !ym   ----> Reverifier lors de la parallelisation des autres schemas
[1632]108
[5101]109  IF (CPPKEY_DEBUGIO) THEN
110    CALL WriteField_u('massem', massem)
111    CALL WriteField_u('wg', wg)
112    CALL WriteField_u('pbarug', pbarug)
113    CALL WriteField_v('pbarvg', pbarvg)
114    CALL WriteField_u('p_tmp', p)
115    CALL WriteField_u('pk_tmp', pk)
116    CALL WriteField_u('teta_tmp', teta)
117    DO iq = 1, nqtot
118      CALL WriteField_u('q_adv' // trim(int2str(iq)), q(:, :, iq))
119    END DO
120  END IF
[1632]121
[5101]122  !  CALL Register_Hallo_v(pbarvg,llm,1,1,1,1,TestRequest)
123  !  CALL SendRequest(TestRequest)
124  !!$OMP BARRIER
125  !  CALL WaitRequest(TestRequest)
126  !$OMP BARRIER
[4058]127
[5101]128  !  WRITE(*,*) 'advtrac 157: appel de vlspltgen_loc'
129  CALL vlspltgen_loc(q, 2., massem, wg, pbarug, pbarvg, dtvr, p, pk, teta)
[1632]130
[5101]131  IF (CPPKEY_DEBUGIO) THEN
132    DO iq = 1, nqtot
133      CALL WriteField_u('q_adv' // trim(int2str(iq)), q(:, :, iq))
134    END DO
135  END IF
136
137  GOTO 1234
138  !-------------------------------------------------------------------------
139  !       Appel des sous programmes d'advection
140  !-------------------------------------------------------------------------
141  DO iq = 1, nqtot
142    !     CALL clock(t_initial)
143    IF(tracers(iq)%parent /= 'air') CYCLE
144    iadv = tracers(iq)%iadv
145    !-----------------------------------------------------------------------
146    SELECT CASE(iadv)
[4056]147      !-----------------------------------------------------------------------
[5101]148    CASE(0); CYCLE
149    !--------------------------------------------------------------------
150    CASE(10)  !--- Schema de Van Leer I MUSCL
151      !--------------------------------------------------------------------
152      !           WRITE(*,*) 'advtrac 239: iq,q(1721,19,:)=',iq,q(1721,19,:)
153      !LF         CALL vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
[1632]154
[5101]155      !--------------------------------------------------------------------
156    CASE(14)  !--- Schema "pseuDO amont" + test sur humidite specifique
157      !--- pour la vapeur d'eau. F. Codron
158      !--------------------------------------------------------------------
159      !           WRITE(*,*) 'advtrac 248: iq,q(1721,19,:)=',iq,q(1721,19,:)
160      CALL abort_gcm("advtrac", "appel a vlspltqs :schema non parallelise", 1)
161      !LF         CALL vlspltqs_p(q(1,1,1),2.,massem,wg,pbarug,pbarvg,dtvr,p,pk,teta )
[1632]162
[5101]163      !--------------------------------------------------------------------
164    CASE(12)  !--- Schema de Frederic Hourdin
165      !--------------------------------------------------------------------
166      CALL abort_gcm("advtrac", "appel a vlspltqs :schema non parallelise", 1)
167      CALL adaptdt(iadv, dtbon, n, pbarug, massem)   ! pas de temps adaptatif
168      IF(n > 1) WRITE(*, *) 'WARNING horizontal dt=', dtbon, 'dtvr=', dtvr, 'n=', n
169      DO indice = 1, n
170        CALL advn(q(1, 1, iq), massem, wg, pbarug, pbarvg, dtbon, 1)
171      END DO
[1632]172
[5101]173      !--------------------------------------------------------------------
174    CASE(13)  !--- Pas de temps adaptatif
175      !--------------------------------------------------------------------
176      CALL abort_gcm("advtrac", "schema non parallelise", 1)
177      CALL adaptdt(iadv, dtbon, n, pbarug, massem)
178      IF(n > 1) WRITE(*, *) 'WARNING horizontal dt=', dtbon, 'dtvr=', dtvr, 'n=', n
179      DO indice = 1, n
180        CALL advn(q(1, 1, iq), massem, wg, pbarug, pbarvg, dtbon, 2)
181      END DO
[1632]182
[5101]183      !--------------------------------------------------------------------
184    CASE(20)  !--- Schema de pente SLOPES
185      !--------------------------------------------------------------------
186      CALL abort_gcm("advtrac", "schema SLOPES non parallelise", 1)
187      CALL pentes_ini (q(1, 1, iq), wg, massem, pbarug, pbarvg, 0)
[1632]188
[5101]189      !--------------------------------------------------------------------
190    CASE(30)  !--- Schema de Prather
191      !--------------------------------------------------------------------
192      CALL abort_gcm("advtrac", "schema prather non parallelise", 1)
193      ! Pas de temps adaptatif
194      CALL adaptdt(iadv, dtbon, n, pbarug, massem)
195      IF(n > 1) WRITE(*, *) 'WARNING horizontal dt=', dtbon, 'dtvr=', dtvr, 'n=', n
196      CALL prather(q(1, 1, iq), wg, massem, pbarug, pbarvg, n, dtbon)
[1632]197
[5101]198      !--------------------------------------------------------------------
199    CASE(11, 16, 17, 18)   !--- Schemas PPM Lin et Rood
200      !--------------------------------------------------------------------
201      CALL abort_gcm("advtrac", "schema PPM non parallelise", 1)
202      ! Test sur le flux horizontal
203      CALL adaptdt(iadv, dtbon, n, pbarug, massem)   ! pas de temps adaptatif
204      IF(n > 1) WRITE(*, *) 'WARNING horizontal dt=', dtbon, 'dtvr=', dtvr, 'n=', n
205      ! Test sur le flux vertical
206      CFLmaxz = 0.
207      DO l = 2, llm
208        DO ij = iip2, ip1jm
209          aaa = wg(ij, l) * dtvr / massem(ij, l)
210          CFLmaxz = max(CFLmaxz, aaa)
211          bbb = -wg(ij, l) * dtvr / massem(ij, l - 1)
212          CFLmaxz = max(CFLmaxz, bbb)
213        END DO
214      END DO
215      IF(CFLmaxz>=1) WRITE(*, *) 'WARNING vertical', 'CFLmaxz=', CFLmaxz
216      !----------------------------------------------------------------
217      !     Ss-prg interface LMDZ.4->PPM3d (ss-prg de Lin)
218      !----------------------------------------------------------------
219      CALL interpre(q(1, 1, iq), qppm(1, 1, iq), wg, fluxwppm, massem, &
220              apppm, bpppm, massebx, masseby, pbarug, pbarvg, &
221              unatppm, vnatppm, psppm)
[1632]222
[5101]223      !----------------------------------------------------------------
224      DO indice = 1, n     !--- VL (version PPM) horiz. et PPM vert.
225        !----------------------------------------------------------------
226        SELECT CASE(iadv)
227          !----------------------------------------------------------
228        CASE(11)
229          !----------------------------------------------------------
230          CALL ppm3d(1, qppm(1, 1, iq), psppm, psppm, unatppm, vnatppm, fluxwppm, dtbon, &
231                  2, 2, 2, 1, iim, jjp1, 2, llm, apppm, bpppm, 0.01, 6400000, fill, dum, 220.)
232          !----------------------------------------------------------
233        CASE(16) !--- Monotonic PPM
234          !----------------------------------------------------------
235          CALL ppm3d(1, qppm(1, 1, iq), psppm, psppm, unatppm, vnatppm, fluxwppm, dtbon, &
236                  3, 3, 3, 1, iim, jjp1, 2, llm, apppm, bpppm, 0.01, 6400000, fill, dum, 220.)
237          !----------------------------------------------------------
238        CASE(17) !--- Semi monotonic PPM
239          !----------------------------------------------------------
240          CALL ppm3d(1, qppm(1, 1, iq), psppm, psppm, unatppm, vnatppm, fluxwppm, dtbon, &
241                  4, 4, 4, 1, iim, jjp1, 2, llm, apppm, bpppm, 0.01, 6400000, fill, dum, 220.)
242          !----------------------------------------------------------
243        CASE(18) !--- Positive Definite PPM
244          !----------------------------------------------------------
245          CALL ppm3d(1, qppm(1, 1, iq), psppm, psppm, unatppm, vnatppm, fluxwppm, dtbon, &
246                  5, 5, 5, 1, iim, jjp1, 2, llm, apppm, bpppm, 0.01, 6400000, fill, dum, 220.)
247        END SELECT
248        !----------------------------------------------------------------
249      END DO
250      !----------------------------------------------------------------
251      !     Ss-prg interface PPM3d-LMDZ.4
252      !----------------------------------------------------------------
253      CALL interpost(q(1, 1, iq), qppm(1, 1, iq))
[4056]254      !----------------------------------------------------------------------
[5101]255    END SELECT
256    !----------------------------------------------------------------------
[1632]257
[5101]258    !----------------------------------------------------------------------
259    ! On impose une seule valeur du traceur au pole Sud j=jjm+1=jjp1 et Nord j=1
260    !----------------------------------------------------------------------
261    !  CALL traceurpole(q(1,1,iq),massem)
[1632]262
[5101]263    !--- Calcul du temps cpu pour un schema donne
264    !  CALL clock(t_final)
265    !ym  tps_cpu=t_final-t_initial
266    !ym  cpuadv(iq)=cpuadv(iq)+tps_cpu
[1632]267
[5101]268  END DO
[1632]269
[5101]270  1234 CONTINUE
271  !$OMP BARRIER
272  IF(planet_type=="earth") THEN
273    ijb = ij_begin
274    ije = ij_end
275    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
276    DO l = 1, llm
277      DO ij = ijb, ije
278        finmasse(ij, l) = p(ij, l) - p(ij, l + 1)
[4056]279      END DO
[5101]280    END DO
281    !$OMP END DO
[1632]282
[5101]283    CALL qminimum_loc(q, nqtot, finmasse)
[1632]284
[5101]285  END IF ! of if (planet_type=="earth")
[1632]286
[4056]287END SUBROUTINE advtrac_loc
[1632]288
Note: See TracBrowser for help on using the repository browser.