source: LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3dpar/advtrac_p.F @ 3867

Last change on this file since 3867 was 1446, checked in by Ehouarn Millour, 14 years ago

Implemented modifications to enable running with only one tracer for planet types different from "earth". Rem: If flag 'planet_type' is set to "earth" (default behaviour) then there must be at least 2 tracers for the dynamics to function properly.

These updates do not induce any changes in model outputs with respect to previous revisions.

EM

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.6 KB
Line 
1!
2! $Id: advtrac_p.F 1446 2010-10-22 09:27:25Z fairhead $
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 ! of IF(iadvtr.EQ.0)
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      if (planet_type=="earth") then
455
456        ijb=ij_begin
457        ije=ij_end
458
459c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
460        DO l = 1, llm
461         DO ij = ijb, ije
462           finmasse(ij,l) =  p(ij,l) - p(ij,l+1)
463         ENDDO
464        ENDDO
465c$OMP END DO
466
467        CALL qminimum_p( q, 2, finmasse )
468
469c------------------------------------------------------------------
470c   on reinitialise a zero les flux de masse cumules
471c---------------------------------------------------
472c          iadvtr=0
473
474c$OMP MASTER
475        call VTe(VTadvection)
476        call stop_timer(timer_vanleer)
477        call VTb(VThallo)
478c$OMP END MASTER
479
480        do j=1,nqtot
481          call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
482     *                             jj_nb_caldyn,0,0,Request_vanleer)
483        enddo
484
485        call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm,
486     *       jj_nb_caldyn,0,0,Request_vanleer)
487
488        call SendRequest(Request_vanleer)
489c$OMP BARRIER
490        call WaitRequest(Request_vanleer)     
491
492c$OMP BARRIER
493c$OMP MASTER
494        call SetDistrib(jj_nb_caldyn)
495        call VTe(VThallo)
496        call resume_timer(timer_caldyn)
497c$OMP END MASTER
498c$OMP BARRIER   
499          iadvtr=0
500      endif ! of if (planet_type=="earth")
501       ENDIF ! if iadvtr.EQ.iapp_tracvl
502
503       RETURN
504       END
505
Note: See TracBrowser for help on using the repository browser.