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

Last change on this file since 1443 was 1441, checked in by emillour, 10 years ago

Updates in common dynamics (seq and ) to keep up with updates
in LMDZ5 (up to LMDZ5 trunk, rev 2250):

  • compilation:
  • added test in grid/dimension/makdim to check that # of longitudes is a multiple of 8
  • dyn3d_common:

Bug correction concerning zoom (cf LMDZ5 rev 2218)

  • coefpoly.F becomes coefpoly_m.F90 (in misc)
  • fxhyp.F => fxhyp_m.F90 , fyhyp.F => fyhyp_m.F90
  • new routines for zoom: invert_zoom_x_m.F90 and principal_cshift_m.F90
  • inigeom.F adapted to new zoom definition routines
  • fluxstokenc.F : got rid of calls to initial0()
  • dyn3d:
  • advtrac.F90 : got rid of calls to initial0()
  • conf_gcm.F90 : cosmetic changes and change in default dzoomx,dzoomy values
  • guide_mod.F90 : followed updates from Earth Model
  • gcm.F is now gcm.F90
  • dyn3dpar:
  • advtrac_p.F90, covcont_p.F90, mod_hallo.F90 : cosmetic changes
  • conf_gcm.F90 : cosmetic and changed in default dzoomx,dzoomy values
  • parallel_lmdz.F90 : updates to keep up with Earth model
  • misc:
  • arth.F90 becomes arth_m.F90
  • wxios.F90 updated wrt Earth model changes
  • nrtype.F90 and coefpoly_m.F90 added
  • ran1.F, sort.F, minmax.F, minmax2.F, juldate.F moved over from dyn3d_common

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 pbaruc(ip1jmp1,llm),pbarvc(ip1jm,llm)
50  REAL massem(ip1jmp1,llm),zdp(ip1jmp1)
51  REAL,SAVE::pbarug(ip1jmp1,llm),pbarvg(ip1jm,llm),wg(ip1jmp1,llm)
52  REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
53  INTEGER iadvtr
54  INTEGER ij,l,iq,iiq
55  REAL zdpmin, zdpmax
56  SAVE iadvtr, massem, pbaruc, pbarvc
57  DATA 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.