source: LMDZ5/trunk/libf/dyn3dpar/advtrac_p.F90 @ 2346

Last change on this file since 2346 was 2239, checked in by Ehouarn Millour, 10 years ago

Reorganizing physics/dynamics interface:

  • what is related to dynamics-physics interface is now in a seperate directory: dynlmdz_phy* for physics in phy*
  • 1d model and related dependencies (including a couple from "dynamics", set up as symbolic links) is now in subdirectory "dyn1d" of phy*.
  • "bibio" directory is now "misc" and should only contain autonomous utilities.
  • "cosp" is now a subdirectory of phylmd.

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