source: LMDZ5/branches/LMDZ5_SPLA/libf/dyn3dpar/advtrac_p.F90 @ 5005

Last change on this file since 5005 was 1987, checked in by Ehouarn Millour, 11 years ago

Add updating pressure, mass and Exner function (ie: all variables which depend on surface pressure) after adding physics tendencies (which include a surface pressure tendency).
Note that this change induces slight changes in GCM results with respect to previous svn version of the code, even if surface pressure tendency is zero (because of recomputation of polar values as an average over polar points on the dynamics grid).
EM
  • 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: 17.0 KB
Line 
1! $Id: advtrac_p.F90 1987 2014-02-24 15:05:47Z dcugnet $
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, ONLY: nqtot, iadv
19  USE control_mod, ONLY: iapp_tracvl, day_step, planet_type
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  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  !-------------------------------------------------------------------
46  !     Ajout PPM
47  !--------------------------------------------------------
48  REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
49  !-------------------------------------------------------------
50  !     Variables locales
51  !-------------------------------------------------------------
52
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)
83
84  ijb_u=ij_begin
85  ije_u=ij_end
86
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
91
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
102
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
114
115  !   selection de la masse instantannee des mailles avant le transport.
116  IF(iadvtr.EQ.0) THEN
117
118     !         CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
119     ijb=ij_begin
120     ije=ij_end
121
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
127
128     !cc         CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
129     !
130  ENDIF ! of IF(iadvtr.EQ.0)
131
132  iadvtr   = iadvtr+1
133
134  !$OMP MASTER
135  iapptrac = iadvtr
136  !$OMP END MASTER
137
138  !   Test pour savoir si on advecte a ce pas de temps
139
140  IF ( iadvtr.EQ.iapp_tracvl ) THEN
141     !$OMP MASTER
142     call suspend_timer(timer_caldyn)
143     !$OMP END MASTER
144
145     ijb=ij_begin
146     ije=ij_end
147
148
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
162        p_tmp(ijb:ije,l)=p(ijb:ije,l)
163     ENDDO
164     !$OMP END DO NOWAIT
165
166     !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
167     DO l=1,llm
168        pk_tmp(ijb:ije,l)=pk(ijb:ije,l)
169        teta_tmp(ijb:ije,l)=teta(ijb:ije,l)
170     ENDDO
171     !$OMP END DO NOWAIT
172
173     !$OMP MASTER
174     call VTb(VTHallo)
175     !$OMP END MASTER
176
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
195
196     call SendRequest(Request_vanleer)
197     !$OMP BARRIER
198     call WaitRequest(Request_vanleer)
199
200
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
209
210     ! ... Flux de masse diaganostiques traceurs
211     ijb=ij_begin
212     ije=ij_end
213     flxw(ijb:ije,1:llm)=wg(ijb:ije,1:llm)/REAL(iapp_tracvl)
214
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
220
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
228
229        !            CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
230        ! ym  ---> pourquoi jjm-1 et non jjm ? a cause du pole ?
231
232        do ij=ijb,ije-iip1+1,iip1
233           zdp(ij)=zdp(ij+iip1-1)
234        enddo
235
236        DO ij = ijb,ije
237           zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
238        ENDDO
239
240
241        !            CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
242        !  ym ---> eventuellement a revoir
243        CALL minmax ( ije-ijb+1, zdp(ijb), zdpmin,zdpmax )
244
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
270     GOTO 1234     
271     !-----------------------------------------------------------
272     !     Appel des sous programmes d'advection
273     !-----------------------------------------------------------
274     do iq=1,nqtot
275        !        call clock(t_initial)
276        if(iadv(iq) == 0) cycle
277        !   ----------------------------------------------------------------
278        !   Schema de Van Leer I MUSCL
279        !   ----------------------------------------------------------------
280        if(iadv(iq).eq.10) THEN
281
282           call vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
283
284           !   ----------------------------------------------------------------
285           !   Schema "pseudo amont" + test sur humidite specifique
286           !    pour la vapeur d'eau. F. Codron
287           !   ----------------------------------------------------------------
288        else if(iadv(iq).eq.14) then
289           !
290           !ym        stop 'advtrac : appel à vlspltqs :schema non parallelise'
291           CALL vlspltqs_p( q(1,1,1), 2., massem, wg , &
292                pbarug,pbarvg,dtvr,p_tmp,pk_tmp,teta_tmp )
293           !   ----------------------------------------------------------------
294           !   Schema de Frederic Hourdin
295           !   ----------------------------------------------------------------
296        else if(iadv(iq).eq.12) then
297           stop 'advtrac : schema non parallelise'
298           !            Pas de temps adaptatif
299           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
300           if (n.GT.1) then
301              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
302                   dtvr,'n=',n
303           endif
304           do indice=1,n
305              call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
306           end do
307        else if(iadv(iq).eq.13) then
308           stop 'advtrac : schema non parallelise'
309           !            Pas de temps adaptatif
310           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
311           if (n.GT.1) then
312              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
313                   dtvr,'n=',n
314           endif
315           do indice=1,n
316              call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
317           end do
318           !   ----------------------------------------------------------------
319           !   Schema de pente SLOPES
320           !   ----------------------------------------------------------------
321        else if (iadv(iq).eq.20) then
322           stop 'advtrac : schema non parallelise'
323
324           call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
325
326           !   ----------------------------------------------------------------
327           !   Schema de Prather
328           !   ----------------------------------------------------------------
329        else if (iadv(iq).eq.30) then
330           stop 'advtrac : schema non parallelise'
331           !            Pas de temps adaptatif
332           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
333           if (n.GT.1) then
334              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
335                   dtvr,'n=',n
336           endif
337           call  prather(q(1,1,iq),wg,massem,pbarug,pbarvg, &
338                n,dtbon)
339           !   ----------------------------------------------------------------
340           !   Schemas PPM Lin et Rood
341           !   ----------------------------------------------------------------
342        else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND. &
343             iadv(iq).LE.18)) then
344
345           stop 'advtrac : schema non parallelise'
346
347           !        Test sur le flux horizontal
348           !        Pas de temps adaptatif
349           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
350           if (n.GT.1) then
351              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
352                   dtvr,'n=',n
353           endif
354           !        Test sur le flux vertical
355           CFLmaxz=0.
356           do l=2,llm
357              do ij=iip2,ip1jm
358                 aaa=wg(ij,l)*dtvr/massem(ij,l)
359                 CFLmaxz=max(CFLmaxz,aaa)
360                 bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
361                 CFLmaxz=max(CFLmaxz,bbb)
362              enddo
363           enddo
364           if (CFLmaxz.GE.1) then
365              write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
366           endif
367
368           !-----------------------------------------------------------
369           !        Ss-prg interface LMDZ.4->PPM3d
370           !-----------------------------------------------------------
371
372           call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, &
373                apppm,bpppm,massebx,masseby,pbarug,pbarvg, &
374                unatppm,vnatppm,psppm)
375
376           do indice=1,n
377              !----------------------------------------------------------------
378              !                         VL (version PPM) horiz. et PPM vert.
379              !----------------------------------------------------------------
380              if (iadv(iq).eq.11) then
381                 !                  Ss-prg PPM3d de Lin
382                 call ppm3d(1,qppm(1,1,iq), &
383                      psppm,psppm, &
384                      unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1, &
385                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
386                      fill,dum,220.)
387
388                 !-------------------------------------------------------------
389                 !                           Monotonic PPM
390                 !-------------------------------------------------------------
391              else if (iadv(iq).eq.16) then
392                 !                  Ss-prg PPM3d de Lin
393                 call ppm3d(1,qppm(1,1,iq), &
394                      psppm,psppm, &
395                      unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1, &
396                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
397                      fill,dum,220.)
398                 !-------------------------------------------------------------
399
400                 !-------------------------------------------------------------
401                 !                           Semi Monotonic PPM
402                 !-------------------------------------------------------------
403              else if (iadv(iq).eq.17) then
404                 !                  Ss-prg PPM3d de Lin
405                 call ppm3d(1,qppm(1,1,iq), &
406                      psppm,psppm, &
407                      unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1, &
408                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
409                      fill,dum,220.)
410                 !-------------------------------------------------------------
411
412                 !-------------------------------------------------------------
413                 !                         Positive Definite PPM
414                 !-------------------------------------------------------------
415              else if (iadv(iq).eq.18) then
416                 !                  Ss-prg PPM3d de Lin
417                 call ppm3d(1,qppm(1,1,iq), &
418                      psppm,psppm, &
419                      unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1, &
420                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
421                      fill,dum,220.)
422                 !-------------------------------------------------------------
423              endif
424           enddo
425           !-----------------------------------------------------------------
426           !               Ss-prg interface PPM3d-LMDZ.4
427           !-----------------------------------------------------------------
428           call interpost(q(1,1,iq),qppm(1,1,iq))
429        endif
430        !----------------------------------------------------------------------
431
432        !-----------------------------------------------------------------
433        ! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
434        ! et Nord j=1
435        !-----------------------------------------------------------------
436
437        !                  call traceurpole(q(1,1,iq),massem)
438
439        ! calcul du temps cpu pour un schema donne
440
441        !                  call clock(t_final)
442        !ym                  tps_cpu=t_final-t_initial
443        !ym                  cpuadv(iq)=cpuadv(iq)+tps_cpu
444
445     end DO
446
4471234 CONTINUE
448     !$OMP BARRIER
449
450     if (planet_type=="earth") then
451
452        ijb=ij_begin
453        ije=ij_end
454
455        !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
456        DO l = 1, llm
457           DO ij = ijb, ije
458              finmasse(ij,l) =  p(ij,l) - p(ij,l+1)
459           ENDDO
460        ENDDO
461        !$OMP END DO
462
463        CALL qminimum_p( q, 2, finmasse )
464     endif ! of if (planet_type=="earth")
465
466        !------------------------------------------------------------------
467        !   on reinitialise a zero les flux de masse cumules
468        !---------------------------------------------------
469        !          iadvtr=0
470
471      !$OMP MASTER
472      call VTe(VTadvection)
473      call stop_timer(timer_vanleer)
474      call VTb(VThallo)
475     !$OMP END MASTER
476
477      do j=1,nqtot
478        call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm, &
479             jj_nb_caldyn,0,0,Request_vanleer)
480      enddo
481
482      call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm, &
483             jj_nb_caldyn,0,0,Request_vanleer)
484
485      call SendRequest(Request_vanleer)
486      !$OMP BARRIER
487      call WaitRequest(Request_vanleer)     
488
489 !$OMP BARRIER
490 !$OMP MASTER
491      call SetDistrib(jj_nb_caldyn)
492      call VTe(VThallo)
493      call resume_timer(timer_caldyn)
494 !$OMP END MASTER
495 !$OMP BARRIER 
496        iadvtr=0
497  ENDIF ! if iadvtr.EQ.iapp_tracvl
498
499END SUBROUTINE advtrac_p
Note: See TracBrowser for help on using the repository browser.