source: LMDZ4/branches/LMDZ4V5.0-LF/libf/dyn3dpar/advtrac_p.F @ 2900

Last change on this file since 2900 was 1299, checked in by Laurent Fairhead, 15 years ago

Nettoyage general pour se rapprocher des normes et éviter des erreurs a la
compilation:

  • tous les FLOAT() sont remplacés par des REAL()
  • tous les STOP dans phylmd sont remplacés par des appels à abort_gcm
  • le common control défini dans le fichier control.h est remplacé par le module control_mod pour éviter des messages sur l'alignement des variables dans les déclarations
  • des $Header$ remplacés par des $Id$ pour svn

Quelques remplacements à faire ont pu m'échapper

General cleanup of the code to try and adhere to norms and to prevent some
compilation errors:

  • all FLOAT() instructions have been replaced by REAL() instructions
  • all STOP instructions in phylmd have been replaced by calls to abort_gcm
  • the common block control defined in the control.h file has been replaced by the control_mod to prevent compilation warnings on the alignement of declared variables
  • $Header$ replaced by $Id$ for svn

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