source: LMDZ5/trunk/libf/dyn3dmem/advtrac_loc.F @ 4095

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

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