source: LMDZ5/branches/IPSLCM5A2.1/libf/dyn3dmem/advtrac_loc.F @ 5017

Last change on this file since 5017 was 2603, checked in by Ehouarn Millour, 8 years ago

Cleanup in the dynamics: turn logic.h into module logic_mod.F90
EM

  • 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: 13.0 KB
Line 
1!
2! $Id$
3!
4c
5c
6#define DEBUG_IO
7#undef DEBUG_IO
8      SUBROUTINE advtrac_loc(pbarug,pbarvg ,wg,
9     *                   p,  massem,q,teta,
10     *                   pk   )
11
12c     Auteur :  F. Hourdin
13c
14c     Modif. P. Le Van     (20/12/97)
15c            F. Codron     (10/99)
16c            D. Le Croller (07/2001)
17c            M.A Filiberti (04/2002)
18c
19      USE parallel_lmdz
20      USE Write_Field_loc
21      USE Write_Field
22      USE Bands
23      USE mod_hallo
24      USE Vampir
25      USE times
26      USE infotrac, ONLY: nqtot, iadv, ok_iso_verif
27      USE control_mod, ONLY: iapp_tracvl, day_step, planet_type
28      USE advtrac_mod, ONLY: finmasse
29      USE comconst_mod, ONLY: dtvr
30      IMPLICIT NONE
31c
32      include "dimensions.h"
33      include "paramet.h"
34      include "comdissip.h"
35      include "comgeom2.h"
36      include "ener.h"
37      include "description.h"
38
39c-------------------------------------------------------------------
40c     Arguments
41c-------------------------------------------------------------------
42c     Ajout PPM
43c--------------------------------------------------------
44      REAL massebx(ijb_u:ije_u,llm),masseby(ijb_v:ije_v,llm)
45c--------------------------------------------------------
46      INTEGER iapptrac
47      REAL pbarug(ijb_u:ije_u,llm),pbarvg(ijb_v:ije_v,llm)
48      REAL wg(ijb_u:ije_u,llm)
49      REAL q(ijb_u:ije_u,llm,nqtot),massem(ijb_u:ije_u,llm)
50      REAL p( ijb_u:ije_u,llmp1 ),teta(ijb_u:ije_u,llm)
51      REAL pk(ijb_u:ije_u,llm)
52
53c-------------------------------------------------------------
54c     Variables locales
55c-------------------------------------------------------------
56
57      REAL zdp(ijb_u:ije_u)
58      REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
59      INTEGER,SAVE :: iadvtr=0
60c$OMP THREADPRIVATE(iadvtr)
61      INTEGER ij,l,iq,iiq
62      REAL zdpmin, zdpmax
63c----------------------------------------------------------
64c     Rajouts pour PPM
65c----------------------------------------------------------
66      INTEGER indice,n
67      REAL dtbon ! Pas de temps adaptatif pour que CFL<1
68      REAL CFLmaxz,aaa,bbb ! CFL maximum
69      REAL psppm(iim,jjb_u:jje_u) ! pression  au sol
70      REAL unatppm(iim,jjb_u:jje_u,llm),vnatppm(iim,jjb_u:jje_u,llm)
71      REAL qppm(iim*jjnb_u,llm,nqtot)
72      REAL fluxwppm(iim,jjb_u:jje_u,llm)
73      REAL apppm(llmp1), bpppm(llmp1)
74      LOGICAL dum,fill
75      DATA fill/.true./
76      DATA dum/.true./
77      integer ijb,ije,ijbu,ijbv,ijeu,ijev,j
78      type(Request),SAVE :: testRequest
79!$OMP THREADPRIVATE(testRequest)
80
81c  test sur l''eventuelle creation de valeurs negatives de la masse
82         ijb=ij_begin
83         ije=ij_end
84         if (pole_nord) ijb=ij_begin+iip1
85         if (pole_sud) ije=ij_end-iip1
86         
87c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
88         DO l=1,llm-1
89            DO ij = ijb+1,ije
90              zdp(ij) =    pbarug(ij-1,l)   - pbarug(ij,l)
91     s                  - pbarvg(ij-iip1,l) + pbarvg(ij,l)
92     s                  +       wg(ij,l+1)  - wg(ij,l)
93            ENDDO
94           
95c            CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
96c ym  ---> pourquoi jjm-1 et non jjm ? a cause du pole ?
97           
98            do ij=ijb,ije-iip1+1,iip1
99              zdp(ij)=zdp(ij+iip1-1)
100            enddo
101           
102            DO ij = ijb,ije
103               zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
104            ENDDO
105
106
107c            CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
108c  ym ---> eventuellement a revoir
109            CALL minmax ( ije-ijb+1, zdp(ijb), zdpmin,zdpmax )
110           
111            IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN
112            PRINT*,'WARNING DP/P l=',l,'  MIN:',zdpmin,
113     s        '   MAX:', zdpmax
114            ENDIF
115
116         ENDDO
117c$OMP END DO NOWAIT
118
119c-------------------------------------------------------------------
120c   Advection proprement dite (Modification Le Croller (07/2001)
121c-------------------------------------------------------------------
122
123c----------------------------------------------------
124c        Calcul des moyennes basées sur la masse
125c----------------------------------------------------
126
127cym      ----> Normalement, inutile pour les schémas classiques
128cym      ----> Revérifier lors de la parallélisation des autres schemas
129   
130cym          call massbar_p(massem,massebx,masseby)         
131
132#ifdef DEBUG_IO   
133          CALL WriteField_u('massem',massem)
134          CALL WriteField_u('wg',wg)
135          CALL WriteField_u('pbarug',pbarug)
136          CALL WriteField_v('pbarvg',pbarvg)
137          CALL WriteField_u('p_tmp',p)
138          CALL WriteField_u('pk_tmp',pk)
139          CALL WriteField_u('teta_tmp',teta)
140          do j=1,nqtot
141            call WriteField_u('q_adv'//trim(int2str(j)),
142     .                q(:,:,j))
143          enddo
144#endif
145
146!         
147!       call Register_Hallo_v(pbarvg,llm,1,1,1,1,TestRequest)
148!
149!       call SendRequest(TestRequest)
150!c$OMP BARRIER
151!       call WaitRequest(TestRequest)
152c$OMP BARRIER
153                 
154          !write(*,*) 'advtrac 157: appel de vlspltgen_loc'
155          call vlspltgen_loc( q,iadv, 2., massem, wg ,
156     *                        pbarug,pbarvg,dtvr,p,
157     *                        pk,teta )
158
159          !write(*,*) 'advtrac 162: apres appel vlspltgen_loc'
160      if (ok_iso_verif) then
161           call check_isotopes(q,ijb_u,ije_u,'advtrac 162')
162      endif !if (ok_iso_verif) then
163
164#ifdef DEBUG_IO     
165          do j=1,nqtot
166            call WriteField_u('q_adv'//trim(int2str(j)),
167     .                q(:,:,j))
168          enddo
169#endif
170         
171          GOTO 1234     
172c-----------------------------------------------------------
173c     Appel des sous programmes d'advection
174c-----------------------------------------------------------
175      do iq=1,nqtot
176c        call clock(t_initial)
177        if(iadv(iq) == 0) cycle
178c   ----------------------------------------------------------------
179c   Schema de Van Leer I MUSCL
180c   ----------------------------------------------------------------
181        if(iadv(iq).eq.10) THEN
182     
183!LF            call vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
184
185c   ----------------------------------------------------------------
186c   Schema "pseudo amont" + test sur humidite specifique
187C    pour la vapeur d'eau. F. Codron
188c   ----------------------------------------------------------------
189        else if(iadv(iq).eq.14) then
190c
191cym           stop 'advtrac : appel à vlspltqs :schema non parallelise'
192!LF           CALL vlspltqs_p( q(1,1,1), 2., massem, wg ,
193!LF     *                 pbarug,pbarvg,dtvr,p,pk,teta )
194c   ----------------------------------------------------------------
195c   Schema de Frederic Hourdin
196c   ----------------------------------------------------------------
197        else if(iadv(iq).eq.12) then
198          stop 'advtrac : schema non parallelise'
199c            Pas de temps adaptatif
200           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
201           if (n.GT.1) then
202           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
203     s             dtvr,'n=',n
204           endif
205           do indice=1,n
206            call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
207           end do
208        else if(iadv(iq).eq.13) then
209          stop 'advtrac : schema non parallelise'
210c            Pas de temps adaptatif
211           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
212           if (n.GT.1) then
213           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
214     s             dtvr,'n=',n
215           endif
216          do indice=1,n
217            call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
218          end do
219c   ----------------------------------------------------------------
220c   Schema de pente SLOPES
221c   ----------------------------------------------------------------
222        else if (iadv(iq).eq.20) then
223          stop 'advtrac : schema non parallelise'
224
225            call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
226
227c   ----------------------------------------------------------------
228c   Schema de Prather
229c   ----------------------------------------------------------------
230        else if (iadv(iq).eq.30) then
231          stop 'advtrac : schema non parallelise'
232c            Pas de temps adaptatif
233           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
234           if (n.GT.1) then
235           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
236     s             dtvr,'n=',n
237           endif
238           call  prather(q(1,1,iq),wg,massem,pbarug,pbarvg,
239     s                     n,dtbon)
240c   ----------------------------------------------------------------
241c   Schemas PPM Lin et Rood
242c   ----------------------------------------------------------------
243       else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND.
244     s                     iadv(iq).LE.18)) then
245
246           stop 'advtrac : schema non parallelise'
247
248c        Test sur le flux horizontal
249c        Pas de temps adaptatif
250         call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
251         if (n.GT.1) then
252           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
253     s             dtvr,'n=',n
254         endif
255c        Test sur le flux vertical
256         CFLmaxz=0.
257         do l=2,llm
258           do ij=iip2,ip1jm
259            aaa=wg(ij,l)*dtvr/massem(ij,l)
260            CFLmaxz=max(CFLmaxz,aaa)
261            bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
262            CFLmaxz=max(CFLmaxz,bbb)
263           enddo
264         enddo
265         if (CFLmaxz.GE.1) then
266            write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
267         endif
268
269c-----------------------------------------------------------
270c        Ss-prg interface LMDZ.4->PPM3d
271c-----------------------------------------------------------
272
273          call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem,
274     s                 apppm,bpppm,massebx,masseby,pbarug,pbarvg,
275     s                 unatppm,vnatppm,psppm)
276
277          do indice=1,n
278c---------------------------------------------------------------------
279c                         VL (version PPM) horiz. et PPM vert.
280c---------------------------------------------------------------------
281                if (iadv(iq).eq.11) then
282c                  Ss-prg PPM3d de Lin
283                  call ppm3d(1,qppm(1,1,iq),
284     s                       psppm,psppm,
285     s                       unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1,
286     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
287     s                       fill,dum,220.)
288
289c----------------------------------------------------------------------
290c                           Monotonic PPM
291c----------------------------------------------------------------------
292               else if (iadv(iq).eq.16) then
293c                  Ss-prg PPM3d de Lin
294                  call ppm3d(1,qppm(1,1,iq),
295     s                       psppm,psppm,
296     s                       unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1,
297     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
298     s                       fill,dum,220.)
299c---------------------------------------------------------------------
300
301c---------------------------------------------------------------------
302c                           Semi Monotonic PPM
303c---------------------------------------------------------------------
304               else if (iadv(iq).eq.17) then
305c                  Ss-prg PPM3d de Lin
306                  call ppm3d(1,qppm(1,1,iq),
307     s                       psppm,psppm,
308     s                       unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1,
309     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
310     s                       fill,dum,220.)
311c---------------------------------------------------------------------
312
313c---------------------------------------------------------------------
314c                         Positive Definite PPM
315c---------------------------------------------------------------------
316                else if (iadv(iq).eq.18) then
317c                  Ss-prg PPM3d de Lin
318                  call ppm3d(1,qppm(1,1,iq),
319     s                       psppm,psppm,
320     s                       unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1,
321     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
322     s                       fill,dum,220.)
323c---------------------------------------------------------------------
324                endif
325            enddo
326c-----------------------------------------------------------------
327c               Ss-prg interface PPM3d-LMDZ.4
328c-----------------------------------------------------------------
329                  call interpost(q(1,1,iq),qppm(1,1,iq))
330            endif
331c----------------------------------------------------------------------
332
333c-----------------------------------------------------------------
334c On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
335c et Nord j=1
336c-----------------------------------------------------------------
337
338c                  call traceurpole(q(1,1,iq),massem)
339
340c calcul du temps cpu pour un schema donne
341
342
343       end DO
344
3451234  CONTINUE
346c$OMP BARRIER
347
348      if (planet_type=="earth") then
349
350      ijb=ij_begin
351      ije=ij_end
352
353c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
354       DO l = 1, llm
355         DO ij = ijb, ije
356           finmasse(ij,l) =  p(ij,l) - p(ij,l+1)
357         ENDDO
358       ENDDO
359c$OMP END DO
360
361        ! CRisi: on passe nqtot et non nq
362       CALL qminimum_loc( q, nqtot, finmasse )
363
364      endif ! of if (planet_type=="earth")
365
366       RETURN
367       END
368
Note: See TracBrowser for help on using the repository browser.