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

Last change on this file since 5120 was 5117, checked in by abarral, 4 months ago

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