source: LMDZ5/trunk/libf/dyn3dpar/advtrac_p.F90 @ 1957

Last change on this file since 1957 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

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