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
Line 
1! $Id: advtrac_p.F90 1549 2011-07-05 08:41:12Z lguez $
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, 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
16  USE mod_hallo
17  USE Vampir
18  USE times
19  USE infotrac, ONLY: nqtot, iadv
20  USE control_mod, ONLY: iapp_tracvl, day_step, planet_type
21  USE comconst_mod, ONLY: dtvr
22  IMPLICIT NONE
23  !
24  include "dimensions.h"
25  include "paramet.h"
26  include "comdissip.h"
27  include "comgeom2.h"
28
29  !-------------------------------------------------------------------
30  !     Arguments
31  !-------------------------------------------------------------------
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  !-------------------------------------------------------------------
42  !     Ajout PPM
43  !--------------------------------------------------------
44  REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
45  !-------------------------------------------------------------
46  !     Variables locales
47  !-------------------------------------------------------------
48
49  REAL,SAVE :: pbaruc(ip1jmp1,llm)=0
50  REAL,SAVE :: pbarvc(ip1jm,llm)=0
51  REAL,SAVE :: massem(ip1jmp1,llm)
52  REAL 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 ij,l,iq,iiq
56  REAL zdpmin, zdpmax
57  INTEGER,SAVE :: iadvtr=0
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)
79
80  ijb_u=ij_begin
81  ije_u=ij_end
82
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
87
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
98
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
110
111  !   selection de la masse instantanee des mailles avant le transport.
112  IF(iadvtr.EQ.0) THEN
113
114     !         CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
115     ijb=ij_begin
116     ije=ij_end
117
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
123
124     !cc         CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
125     !
126  ENDIF ! of IF(iadvtr.EQ.0)
127
128  iadvtr   = iadvtr+1
129
130  !$OMP MASTER
131  iapptrac = iadvtr
132  !$OMP END MASTER
133
134  !   Test pour savoir si on advecte a ce pas de temps
135
136  IF ( iadvtr.EQ.iapp_tracvl ) THEN
137     !$OMP MASTER
138     call suspend_timer(timer_caldyn)
139     !$OMP END MASTER
140
141     ijb=ij_begin
142     ije=ij_end
143
144
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
158        p_tmp(ijb:ije,l)=p(ijb:ije,l)
159     ENDDO
160     !$OMP END DO NOWAIT
161
162     !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
163     DO l=1,llm
164        pk_tmp(ijb:ije,l)=pk(ijb:ije,l)
165        teta_tmp(ijb:ije,l)=teta(ijb:ije,l)
166     ENDDO
167     !$OMP END DO NOWAIT
168
169     !$OMP MASTER
170     call VTb(VTHallo)
171     !$OMP END MASTER
172
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
191
192     call SendRequest(Request_vanleer)
193     !$OMP BARRIER
194     call WaitRequest(Request_vanleer)
195
196
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
205
206     ! ... Flux de masse diagnostiques traceurs
207     ijb=ij_begin
208     ije=ij_end
209     flxw(ijb:ije,1:llm)=wg(ijb:ije,1:llm)/REAL(iapp_tracvl)
210
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
216
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
224
225        !            CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
226        ! ym  ---> pourquoi jjm-1 et non jjm ? a cause du pole ?
227
228        do ij=ijb,ije-iip1+1,iip1
229           zdp(ij)=zdp(ij+iip1-1)
230        enddo
231
232        DO ij = ijb,ije
233           zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
234        ENDDO
235
236
237        !            CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
238        !  ym ---> eventuellement a revoir
239        CALL minmax ( ije-ijb+1, zdp(ijb), zdpmin,zdpmax )
240
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
266!!! ATTENTION !!!! TOUT CE QUI EST ENTRE ICI ET 1234 EST OBSOLETE !!!!!!!
267     GOTO 1234
268!!! ATTENTION !!!!
269
270     !-----------------------------------------------------------
271     !     Appel des sous programmes d'advection
272     !-----------------------------------------------------------
273     do iq=1,nqtot
274        !        call clock(t_initial)
275        if(iadv(iq) == 0) cycle
276        !   ----------------------------------------------------------------
277        !   Schema de Van Leer I MUSCL
278        !   ----------------------------------------------------------------
279        if(iadv(iq).eq.10) THEN
280
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           !   ----------------------------------------------------------------
287        else if(iadv(iq).eq.14) then
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           !   ----------------------------------------------------------------
295        else if(iadv(iq).eq.12) then
296           stop 'advtrac : schema non parallelise'
297           !            Pas de temps adaptatif
298           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
299           if (n.GT.1) then
300              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
301                   dtvr,'n=',n
302           endif
303           do indice=1,n
304              call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
305           end do
306        else if(iadv(iq).eq.13) then
307           stop 'advtrac : schema non parallelise'
308           !            Pas de temps adaptatif
309           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
310           if (n.GT.1) then
311              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
312                   dtvr,'n=',n
313           endif
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           !   ----------------------------------------------------------------
320        else if (iadv(iq).eq.20) then
321           stop 'advtrac : schema non parallelise'
322
323           call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
324
325           !   ----------------------------------------------------------------
326           !   Schema de Prather
327           !   ----------------------------------------------------------------
328        else if (iadv(iq).eq.30) then
329           stop 'advtrac : schema non parallelise'
330           !            Pas de temps adaptatif
331           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
332           if (n.GT.1) then
333              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
334                   dtvr,'n=',n
335           endif
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
343
344           stop 'advtrac : schema non parallelise'
345
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
362           enddo
363           if (CFLmaxz.GE.1) then
364              write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
365           endif
366
367           !-----------------------------------------------------------
368           !        Ss-prg interface LMDZ.4->PPM3d
369           !-----------------------------------------------------------
370
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)
374
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.)
386
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                 !-------------------------------------------------------------
398
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                 !-------------------------------------------------------------
410
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        !----------------------------------------------------------------------
430
431        !-----------------------------------------------------------------
432        ! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
433        ! et Nord j=1
434        !-----------------------------------------------------------------
435
436        !                  call traceurpole(q(1,1,iq),massem)
437
438        ! calcul du temps cpu pour un schema donne
439
440        !                  call clock(t_final)
441        !ym                  tps_cpu=t_final-t_initial
442        !ym                  cpuadv(iq)=cpuadv(iq)+tps_cpu
443
444     end DO
445
446!!! ATTENTION !!!!
4471234 CONTINUE
448!!! ATTENTION !!!! LE CODE REPREND ICI !!!!!!!!
449
450     !$OMP BARRIER
451
452     if (planet_type=="earth") then
453
454        ijb=ij_begin
455        ije=ij_end
456
457        !$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
463        !$OMP END DO
464
465        CALL qminimum_p( q, 2, finmasse )
466
467     endif ! of if (planet_type=="earth")
468
469        !------------------------------------------------------------------
470        !   on reinitialise a zero les flux de masse cumules
471        !---------------------------------------------------
472        !          iadvtr=0
473
474        !$OMP MASTER
475        call VTe(VTadvection)
476        call stop_timer(timer_vanleer)
477        call VTb(VThallo)
478        !$OMP END MASTER
479
480
481        do j=1,nqtot
482           call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm, &
483                jj_nb_caldyn,0,0,Request_vanleer)
484        enddo
485
486        call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm, &
487             jj_nb_caldyn,0,0,Request_vanleer)
488
489        call SendRequest(Request_vanleer)
490        !$OMP BARRIER
491        call WaitRequest(Request_vanleer)     
492
493        !$OMP BARRIER
494        !$OMP MASTER
495        call SetDistrib(jj_nb_caldyn)
496        call VTe(VThallo)
497        call resume_timer(timer_caldyn)
498 !$OMP END MASTER
499 !$OMP BARRIER 
500        iadvtr=0
501  ENDIF ! if iadvtr.EQ.iapp_tracvl
502
503END SUBROUTINE advtrac_p
Note: See TracBrowser for help on using the repository browser.