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

Last change on this file since 2014 was 1959, checked in by emillour, 7 years ago

Common (parallel) dynamics:
Add some initializations of global arrays. This is not a bug fix, strictly
speaking, as results are unchanged. But in some cases computations extend
into halos and these cases are identified as errors with recent versions
of ifort (2018) with the -init=arrays,snan debugging option.
EM

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