source: trunk/LMDZ.COMMON/libf/dyn3dpar/advtrac_p.F90 @ 1231

Last change on this file since 1231 was 1189, checked in by emillour, 11 years ago

Common dynamics: a couple of bug fixes:

  • Correctly account for the change in pressure, mass, etc. after modifying surface pressure following a call to the physics.
  • Corrected tracer advection, which is computed using values at the beginning of the time step, so it is done at Matsuno forward step.

EM

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