source: LMDZ4/branches/V3_test/libf/dyn3dpar/advtrac_p.F @ 718

Last change on this file since 718 was 709, checked in by Laurent Fairhead, 18 years ago

Nouvelles versions de la dynamique YM
LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 17.0 KB
RevLine 
[630]1!
2! $Header$
3!
4c
5c
[709]6#ifdef INCA
[630]7      SUBROUTINE advtrac_p(pbaru,pbarv ,
8     *                   p,  masse,q,iapptrac,teta,
9     *                  flxw,
10     *                  pk,
11     *                  mmt_adj,
12     *                  hadv_flg)
13#else
14      SUBROUTINE advtrac_p(pbaru,pbarv ,
15     *                   p,  masse,q,iapptrac,teta,
16     *                  pk)
17#endif
18
19c     Auteur :  F. Hourdin
20c
21c     Modif. P. Le Van     (20/12/97)
22c            F. Codron     (10/99)
23c            D. Le Croller (07/2001)
24c            M.A Filiberti (04/2002)
25c
26      USE parallel
[709]27      USE Write_Field_p
[630]28      USE Bands
29      USE mod_hallo
30      USE Vampir
31      USE times
32      IMPLICIT NONE
33c
34#include "dimensions.h"
35#include "paramet.h"
36#include "comconst.h"
37#include "comvert.h"
38#include "comdissip.h"
39#include "comgeom2.h"
40#include "logic.h"
41#include "temps.h"
42#include "control.h"
43#include "ener.h"
44#include "description.h"
45#include "advtrac.h"
46
47c-------------------------------------------------------------------
48c     Arguments
49c-------------------------------------------------------------------
50c     Ajout PPM
51c--------------------------------------------------------
52      REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
53c--------------------------------------------------------
54      INTEGER iapptrac
55      REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
56      REAL q(ip1jmp1,llm,nqmx),masse(ip1jmp1,llm)
57      REAL p( ip1jmp1,llmp1 ),teta(ip1jmp1,llm)
58      REAL pk(ip1jmp1,llm)
[709]59#ifdef INCA
[630]60cym      INTEGER            :: hadv_flg(nq)
61         INTEGER            :: hadv_flg(nqmx)
62cym      REAL               :: mmt_adj(ip1jmp1,llm)
63       REAL               :: mmt_adj(ip1jmp1,llm,1)
64       REAL               :: flxw(ip1jmp1,llm)
65#endif
66
67c-------------------------------------------------------------
68c     Variables locales
69c-------------------------------------------------------------
70
71      REAL pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
72      REAL massem(ip1jmp1,llm),zdp(ip1jmp1)
[709]73      REAL,SAVE::pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm)
[630]74      REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
75      real cpuadv(nqmx)
76      common/cpuadv/cpuadv
77
78      INTEGER iadvtr
79      INTEGER ij,l,iq,iiq
80      REAL zdpmin, zdpmax
81      SAVE iadvtr, massem, pbaruc, pbarvc
82      DATA iadvtr/0/
[709]83c$OMP THREADPRIVATE(iadvtr)
[630]84c----------------------------------------------------------
85c     Rajouts pour PPM
86c----------------------------------------------------------
87      INTEGER indice,n
88      REAL dtbon ! Pas de temps adaptatif pour que CFL<1
89      REAL CFLmaxz,aaa,bbb ! CFL maximum
90      REAL psppm(iim,jjp1) ! pression  au sol
91      REAL unatppm(iim,jjp1,llm),vnatppm(iim,jjp1,llm)
92      REAL qppm(iim*jjp1,llm,nqmx)
93      REAL fluxwppm(iim,jjp1,llm)
94      REAL apppm(llmp1), bpppm(llmp1)
95      LOGICAL dum,fill
96      DATA fill/.true./
97      DATA dum/.true./
98      REAL finmasse(ip1jmp1,llm)
99      integer ijb,ije,ijb_u,ijb_v,ije_u,ije_v,j
100      type(Request) :: Request_vanleer
[709]101      REAL,SAVE :: p_tmp( ip1jmp1,llmp1 )
102      REAL,SAVE :: teta_tmp(ip1jmp1,llm)
103      REAL,SAVE :: pk_tmp(ip1jmp1,llm)
[630]104     
105      ijb_u=ij_begin
106      ije_u=ij_end
107     
108      ijb_v=ij_begin-iip1
109      ije_v=ij_end
110      if (pole_nord) ijb_v=ij_begin
111      if (pole_sud)  ije_v=ij_end-iip1
112
113      IF(iadvtr.EQ.0) THEN
114c         CALL initial0(ijp1llm,pbaruc)
115c         CALL initial0(ijmllm,pbarvc)
[709]116c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
117        DO l=1,llm   
118          pbaruc(ijb_u:ije_u,l)=0.
119          pbarvc(ijb_v:ije_v,l)=0.
120        ENDDO
121c$OMP END DO NOWAIT 
[630]122      ENDIF
123
124c   accumulation des flux de masse horizontaux
[709]125c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[630]126      DO l=1,llm
127         DO ij = ijb_u,ije_u
128            pbaruc(ij,l) = pbaruc(ij,l) + pbaru(ij,l)
129         ENDDO
130         DO ij = ijb_v,ije_v
131            pbarvc(ij,l) = pbarvc(ij,l) + pbarv(ij,l)
132         ENDDO
133      ENDDO
[709]134c$OMP END DO NOWAIT
[630]135
136c   selection de la masse instantannee des mailles avant le transport.
137      IF(iadvtr.EQ.0) THEN
138
139c         CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
140          ijb=ij_begin
141          ije=ij_end
142
[709]143c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
144       DO l=1,llm
145          massem(ijb:ije,l)=masse(ijb:ije,l)
146       ENDDO
147c$OMP END DO NOWAIT
148
[630]149ccc         CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
150c
151      ENDIF
152
153      iadvtr   = iadvtr+1
[709]154
155c$OMP MASTER
[630]156      iapptrac = iadvtr
[709]157c$OMP END MASTER
[630]158
[709]159c   Test pour savoir si on advecte a ce pas de temps
[630]160
161      IF ( iadvtr.EQ.iapp_tracvl ) THEN
[709]162c$OMP MASTER
[630]163        call suspend_timer(timer_caldyn)
[709]164c$OMP END MASTER
[630]165     
166      ijb=ij_begin
167      ije=ij_end
168     
169
170cc   ..  Modif P.Le Van  ( 20/12/97 )  ....
171cc
172
173c   traitement des flux de masse avant advection.
174c     1. calcul de w
175c     2. groupement des mailles pres du pole.
176
[709]177c$OMP BARRIER
[630]178        CALL groupe_p( massem, pbaruc,pbarvc, pbarug,pbarvg,wg )
[709]179c$OMP BARRIER
[630]180
[709]181c$OMP BARRIER
182c$OMP MASTER     
183      p_tmp(ijb:ije,1:llmp1)=p(ijb:ije,1:llmp1)
184      pk_tmp(ijb:ije,1:llm)=pk(ijb:ije,1:llm)
185      teta_tmp(ijb:ije,1:llm)=teta(ijb:ije,1:llm)
[630]186      call VTb(VTHallo)
187      call Register_SwapFieldHallo(pbarug,pbarug,ip1jmp1,llm,
188     *                             jj_Nb_vanleer,0,0,Request_vanleer)
189      call Register_SwapFieldHallo(pbarvg,pbarvg,ip1jm,llm,
190     *                             jj_Nb_vanleer,1,0,Request_vanleer)
191      call Register_SwapFieldHallo(massem,massem,ip1jmp1,llm,
192     *                             jj_Nb_vanleer,0,0,Request_vanleer)
193      call Register_SwapFieldHallo(wg,wg,ip1jmp1,llm,
194     *                             jj_Nb_vanleer,0,0,Request_vanleer)
195      call Register_SwapFieldHallo(teta_tmp,teta_tmp,ip1jmp1,llm,
196     *                             jj_Nb_vanleer,1,1,Request_vanleer)
197      call Register_SwapFieldHallo(p_tmp,p_tmp,ip1jmp1,llmp1,
198     *                             jj_Nb_vanleer,1,1,Request_vanleer)
199      call Register_SwapFieldHallo(pk_tmp,pk_tmp,ip1jmp1,llm,
200     *                             jj_Nb_vanleer,1,1,Request_vanleer)
201      do j=1,nqmx
202        call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
203     *                             jj_nb_vanleer,0,0,Request_vanleer)
204      enddo
205     
206      call SetDistrib(jj_nb_vanleer)
207      call SendRequest(Request_vanleer)
208      call WaitRequest(Request_vanleer)
209
210      call VTe(VTHallo)
211     
212      call VTb(VTadvection)
213      call start_timer(timer_vanleer)
[709]214c$OMP END MASTER
215c$OMP BARRIER
[630]216     
[709]217#ifdef INCA
[630]218      ! ... Flux de masse diaganostiques traceurs
219c      flxw = wg / FLOAT(iapp_tracvl)
220         ijb=ij_begin
221         ije=ij_end
222         flxw(ijb:ije,1:llm)=wg(ijb:ije,1:llm)/FLOAT(iapp_tracvl)
223#endif
224
225c  test sur l'eventuelle creation de valeurs negatives de la masse
226         ijb=ij_begin
227         ije=ij_end
228         if (pole_nord) ijb=ij_begin+iip1
229         if (pole_sud) ije=ij_end-iip1
230         
[709]231c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
[630]232         DO l=1,llm-1
233            DO ij = ijb+1,ije
234              zdp(ij) =    pbarug(ij-1,l)   - pbarug(ij,l)
235     s                  - pbarvg(ij-iip1,l) + pbarvg(ij,l)
236     s                  +       wg(ij,l+1)  - wg(ij,l)
237            ENDDO
238           
239c            CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
240c ym  ---> pourquoi jjm-1 et non jjm ? a cause du pole ?
241           
242            do ij=ijb,ije-iip1+1,iip1
243              zdp(ij)=zdp(ij+iip1-1)
244            enddo
245           
246            DO ij = ijb,ije
247               zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
248            ENDDO
249
250
251c            CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
252c  ym ---> eventuellement a revoir
253            CALL minmax ( ije-ijb+1, zdp(ijb), zdpmin,zdpmax )
254           
255            IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN
256            PRINT*,'WARNING DP/P l=',l,'  MIN:',zdpmin,
257     s        '   MAX:', zdpmax
258            ENDIF
259
260         ENDDO
[709]261c$OMP END DO NOWAIT
[630]262
263c-------------------------------------------------------------------
264c   Advection proprement dite (Modification Le Croller (07/2001)
265c-------------------------------------------------------------------
266
267c----------------------------------------------------
268c        Calcul des moyennes basées sur la masse
269c----------------------------------------------------
270
271cym      ----> Normalement, inutile pour les schémas classiques
272cym      ----> Revérifier lors de la parallélisation des autres schemas
273   
[709]274cym          call massbar_p(massem,massebx,masseby)         
[630]275
[709]276          call vlspltgen_p( q,iadv, 2., massem, wg ,
277     *                    pbarug,pbarvg,dtvr,p_tmp,pk_tmp,teta_tmp )
278
279         
280          GOTO 1234     
[630]281c-----------------------------------------------------------
282c     Appel des sous programmes d'advection
283c-----------------------------------------------------------
284      do iq=1,nqmx
285c        call clock(t_initial)
286        if(iadv(iq) == 0) cycle
287c   ----------------------------------------------------------------
288c   Schema de Van Leer I MUSCL
289c   ----------------------------------------------------------------
290        if(iadv(iq).eq.10) THEN
291     
292            call vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
293
294c   ----------------------------------------------------------------
295c   Schema "pseudo amont" + test sur humidite specifique
296C    pour la vapeur d'eau. F. Codron
297c   ----------------------------------------------------------------
298        else if(iadv(iq).eq.14) then
299c
300cym           stop 'advtrac : appel à vlspltqs :schema non parallelise'
301           CALL vlspltqs_p( q(1,1,1), 2., massem, wg ,
302     *                 pbarug,pbarvg,dtvr,p_tmp,pk_tmp,teta_tmp )
303c   ----------------------------------------------------------------
304c   Schema de Frederic Hourdin
305c   ----------------------------------------------------------------
306        else if(iadv(iq).eq.12) then
307          stop 'advtrac : schema non parallelise'
308c            Pas de temps adaptatif
309           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
310           if (n.GT.1) then
311           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
312     s             dtvr,'n=',n
313           endif
314           do indice=1,n
315            call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
316           end do
317        else if(iadv(iq).eq.13) then
318          stop 'advtrac : schema non parallelise'
319c            Pas de temps adaptatif
320           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
321           if (n.GT.1) then
322           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
323     s             dtvr,'n=',n
324           endif
325          do indice=1,n
326            call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
327          end do
328c   ----------------------------------------------------------------
329c   Schema de pente SLOPES
330c   ----------------------------------------------------------------
331        else if (iadv(iq).eq.20) then
332          stop 'advtrac : schema non parallelise'
333
334            call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
[709]335#ifdef INCA
[630]336       do iiq = iq+1, iq+3
337         q(:,:,iiq)=q(:,:,iiq)*mmt_adj(:,:,1)
338       enddo
339#endif
340
341c   ----------------------------------------------------------------
342c   Schema de Prather
343c   ----------------------------------------------------------------
344        else if (iadv(iq).eq.30) then
345          stop 'advtrac : schema non parallelise'
346c            Pas de temps adaptatif
347           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
348           if (n.GT.1) then
349           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
350     s             dtvr,'n=',n
351           endif
352           call  prather(q(1,1,iq),wg,massem,pbarug,pbarvg,
353     s                     n,dtbon)
[709]354#ifdef INCA
[630]355       do iiq = iq+1, iq+9
356         q(:,:,iiq)=q(:,:,iiq)*mmt_adj(:,:,1)
357       enddo
358#endif
359c   ----------------------------------------------------------------
360c   Schemas PPM Lin et Rood
361c   ----------------------------------------------------------------
362       else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND.
363     s                     iadv(iq).LE.18)) then
364
365           stop 'advtrac : schema non parallelise'
366
367c        Test sur le flux horizontal
368c        Pas de temps adaptatif
369         call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
370         if (n.GT.1) then
371           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
372     s             dtvr,'n=',n
373         endif
374c        Test sur le flux vertical
375         CFLmaxz=0.
376         do l=2,llm
377           do ij=iip2,ip1jm
378            aaa=wg(ij,l)*dtvr/massem(ij,l)
379            CFLmaxz=max(CFLmaxz,aaa)
380            bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
381            CFLmaxz=max(CFLmaxz,bbb)
382           enddo
383         enddo
384         if (CFLmaxz.GE.1) then
385            write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
386         endif
387
388c-----------------------------------------------------------
389c        Ss-prg interface LMDZ.4->PPM3d
390c-----------------------------------------------------------
391
392          call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem,
393     s                 apppm,bpppm,massebx,masseby,pbarug,pbarvg,
394     s                 unatppm,vnatppm,psppm)
395
396          do indice=1,n
397c---------------------------------------------------------------------
398c                         VL (version PPM) horiz. et PPM vert.
399c---------------------------------------------------------------------
400                if (iadv(iq).eq.11) then
401c                  Ss-prg PPM3d de Lin
402                  call ppm3d(1,qppm(1,1,iq),
403     s                       psppm,psppm,
404     s                       unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1,
405     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
406     s                       fill,dum,220.)
407
408c----------------------------------------------------------------------
409c                           Monotonic PPM
410c----------------------------------------------------------------------
411               else if (iadv(iq).eq.16) then
412c                  Ss-prg PPM3d de Lin
413                  call ppm3d(1,qppm(1,1,iq),
414     s                       psppm,psppm,
415     s                       unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1,
416     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
417     s                       fill,dum,220.)
418c---------------------------------------------------------------------
419
420c---------------------------------------------------------------------
421c                           Semi Monotonic PPM
422c---------------------------------------------------------------------
423               else if (iadv(iq).eq.17) then
424c                  Ss-prg PPM3d de Lin
425                  call ppm3d(1,qppm(1,1,iq),
426     s                       psppm,psppm,
427     s                       unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1,
428     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
429     s                       fill,dum,220.)
430c---------------------------------------------------------------------
431
432c---------------------------------------------------------------------
433c                         Positive Definite PPM
434c---------------------------------------------------------------------
435                else if (iadv(iq).eq.18) then
436c                  Ss-prg PPM3d de Lin
437                  call ppm3d(1,qppm(1,1,iq),
438     s                       psppm,psppm,
439     s                       unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1,
440     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
441     s                       fill,dum,220.)
442c---------------------------------------------------------------------
443                endif
444            enddo
445c-----------------------------------------------------------------
446c               Ss-prg interface PPM3d-LMDZ.4
447c-----------------------------------------------------------------
448                  call interpost(q(1,1,iq),qppm(1,1,iq))
449            endif
450c----------------------------------------------------------------------
451
452c-----------------------------------------------------------------
453c On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
454c et Nord j=1
455c-----------------------------------------------------------------
456
457c                  call traceurpole(q(1,1,iq),massem)
458
459c calcul du temps cpu pour un schema donne
460
461c                  call clock(t_final)
462cym                  tps_cpu=t_final-t_initial
463cym                  cpuadv(iq)=cpuadv(iq)+tps_cpu
464
465       end DO
466
[709]4671234  CONTINUE
468c$OMP BARRIER
469c$OMP MASTER
[630]470      ijb=ij_begin
471      ije=ij_end
472     
473       DO l = 1, llm
474         DO ij = ijb, ije
475           finmasse(ij,l) =  p(ij,l) - p(ij,l+1)
476         ENDDO
477       ENDDO
478
[709]479
[630]480       CALL qminimum_p( q, 2, finmasse )
481
482c------------------------------------------------------------------
483c   on reinitialise a zero les flux de masse cumules
484c---------------------------------------------------
[709]485c          iadvtr=0
[630]486        call VTe(VTadvection)
487        call stop_timer(timer_vanleer)
488
489        call VTb(VThallo)
490        do j=1,nqmx
491          call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
492     *                             jj_nb_caldyn,0,0,Request_vanleer)
493        enddo
494
[709]495#ifdef INCA
[630]496       call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm,
497     *                             jj_nb_caldyn,0,0,Request_vanleer)
498#endif     
499        call SetDistrib(jj_nb_caldyn)
500        call SendRequest(Request_vanleer)
501        call WaitRequest(Request_vanleer)     
502       
503        call VTe(VThallo)
504        call resume_timer(timer_caldyn)
[709]505c$OMP END MASTER
506c$OMP BARRIER   
507          iadvtr=0
[630]508       ENDIF ! if iadvtr.EQ.iapp_tracvl
509
510       RETURN
511       END
512
Note: See TracBrowser for help on using the repository browser.