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

Last change on this file since 5172 was 5159, checked in by abarral, 6 months ago

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