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

Last change on this file since 2375 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
RevLine 
[1403]1! $Id: advtrac_p.F90 2239 2015-03-23 07:27:30Z emillour $
[630]2
[1549]3SUBROUTINE advtrac_p(pbaru,pbarv , p,  masse,q,iapptrac,teta, flxw, pk)
[630]4
[1549]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  !
[2239]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
[1549]16  USE mod_hallo
17  USE Vampir
18  USE times
[1987]19  USE infotrac, ONLY: nqtot, iadv
20  USE control_mod, ONLY: iapp_tracvl, day_step, planet_type
[1549]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"
[630]33
[1549]34  !-------------------------------------------------------------------
35  !     Arguments
36  !-------------------------------------------------------------------
[1987]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  !-------------------------------------------------------------------
[1549]47  !     Ajout PPM
48  !--------------------------------------------------------
49  REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm)
50  !-------------------------------------------------------------
51  !     Variables locales
52  !-------------------------------------------------------------
[1146]53
[1549]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)
[630]84
[1549]85  ijb_u=ij_begin
86  ije_u=ij_end
[630]87
[1549]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
[630]92
[1549]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
[630]103
[1549]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
[630]115
[1549]116  !   selection de la masse instantannee des mailles avant le transport.
117  IF(iadvtr.EQ.0) THEN
[764]118
[1549]119     !         CALL SCOPY(ip1jmp1*llm,masse,1,massem,1)
120     ijb=ij_begin
121     ije=ij_end
[630]122
[1549]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
[764]128
[1549]129     !cc         CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
130     !
131  ENDIF ! of IF(iadvtr.EQ.0)
[630]132
[1549]133  iadvtr   = iadvtr+1
[630]134
[1549]135  !$OMP MASTER
136  iapptrac = iadvtr
137  !$OMP END MASTER
[630]138
[1549]139  !   Test pour savoir si on advecte a ce pas de temps
[630]140
[1549]141  IF ( iadvtr.EQ.iapp_tracvl ) THEN
142     !$OMP MASTER
143     call suspend_timer(timer_caldyn)
144     !$OMP END MASTER
[630]145
[1549]146     ijb=ij_begin
147     ije=ij_end
[985]148
[630]149
[1549]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
[985]163        p_tmp(ijb:ije,l)=p(ijb:ije,l)
[1549]164     ENDDO
165     !$OMP END DO NOWAIT
166
167     !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
168     DO l=1,llm
[985]169        pk_tmp(ijb:ije,l)=pk(ijb:ije,l)
170        teta_tmp(ijb:ije,l)=teta(ijb:ije,l)
[1549]171     ENDDO
172     !$OMP END DO NOWAIT
[985]173
[1549]174     !$OMP MASTER
175     call VTb(VTHallo)
176     !$OMP END MASTER
[985]177
[1549]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
[985]196
[1549]197     call SendRequest(Request_vanleer)
198     !$OMP BARRIER
199     call WaitRequest(Request_vanleer)
[630]200
[985]201
[1549]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
[630]210
[1549]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)
[630]215
[1549]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
[630]221
[1549]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
[630]229
[1549]230        !            CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
231        ! ym  ---> pourquoi jjm-1 et non jjm ? a cause du pole ?
[630]232
[1549]233        do ij=ijb,ije-iip1+1,iip1
234           zdp(ij)=zdp(ij+iip1-1)
235        enddo
[630]236
[1549]237        DO ij = ijb,ije
238           zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
239        ENDDO
[630]240
241
[1549]242        !            CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
243        !  ym ---> eventuellement a revoir
244        CALL minmax ( ije-ijb+1, zdp(ijb), zdpmin,zdpmax )
[764]245
[1549]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)
[630]277        if(iadv(iq) == 0) cycle
[1549]278        !   ----------------------------------------------------------------
279        !   Schema de Van Leer I MUSCL
280        !   ----------------------------------------------------------------
[630]281        if(iadv(iq).eq.10) THEN
282
[1549]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           !   ----------------------------------------------------------------
[630]289        else if(iadv(iq).eq.14) then
[1549]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           !   ----------------------------------------------------------------
[630]297        else if(iadv(iq).eq.12) then
[1549]298           stop 'advtrac : schema non parallelise'
299           !            Pas de temps adaptatif
[630]300           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
301           if (n.GT.1) then
[1549]302              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
303                   dtvr,'n=',n
[630]304           endif
305           do indice=1,n
[1549]306              call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
[630]307           end do
308        else if(iadv(iq).eq.13) then
[1549]309           stop 'advtrac : schema non parallelise'
310           !            Pas de temps adaptatif
[630]311           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
312           if (n.GT.1) then
[1549]313              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
314                   dtvr,'n=',n
[630]315           endif
[1549]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           !   ----------------------------------------------------------------
[630]322        else if (iadv(iq).eq.20) then
[1549]323           stop 'advtrac : schema non parallelise'
[630]324
[1549]325           call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
[630]326
[1549]327           !   ----------------------------------------------------------------
328           !   Schema de Prather
329           !   ----------------------------------------------------------------
[630]330        else if (iadv(iq).eq.30) then
[1549]331           stop 'advtrac : schema non parallelise'
332           !            Pas de temps adaptatif
[630]333           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
334           if (n.GT.1) then
[1549]335              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
336                   dtvr,'n=',n
[630]337           endif
[1549]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
[630]345
346           stop 'advtrac : schema non parallelise'
347
[1549]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
[630]364           enddo
[1549]365           if (CFLmaxz.GE.1) then
366              write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
367           endif
[630]368
[1549]369           !-----------------------------------------------------------
370           !        Ss-prg interface LMDZ.4->PPM3d
371           !-----------------------------------------------------------
[630]372
[1549]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)
[630]376
[1549]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.)
[630]388
[1549]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                 !-------------------------------------------------------------
[630]400
[1549]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                 !-------------------------------------------------------------
[630]412
[1549]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        !----------------------------------------------------------------------
[630]432
[1549]433        !-----------------------------------------------------------------
434        ! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
435        ! et Nord j=1
436        !-----------------------------------------------------------------
[630]437
[1549]438        !                  call traceurpole(q(1,1,iq),massem)
[630]439
[1549]440        ! calcul du temps cpu pour un schema donne
[630]441
[1549]442        !                  call clock(t_final)
443        !ym                  tps_cpu=t_final-t_initial
444        !ym                  cpuadv(iq)=cpuadv(iq)+tps_cpu
[630]445
[1549]446     end DO
[630]447
[1549]4481234 CONTINUE
449     !$OMP BARRIER
[985]450
[1549]451     if (planet_type=="earth") then
[985]452
[1454]453        ijb=ij_begin
454        ije=ij_end
455
[1549]456        !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
[1454]457        DO l = 1, llm
[1549]458           DO ij = ijb, ije
459              finmasse(ij,l) =  p(ij,l) - p(ij,l+1)
460           ENDDO
[1454]461        ENDDO
[1549]462        !$OMP END DO
[630]463
[1454]464        CALL qminimum_p( q, 2, finmasse )
[1844]465     endif ! of if (planet_type=="earth")
[630]466
[1549]467        !------------------------------------------------------------------
468        !   on reinitialise a zero les flux de masse cumules
469        !---------------------------------------------------
470        !          iadvtr=0
[985]471
[1844]472      !$OMP MASTER
473      call VTe(VTadvection)
474      call stop_timer(timer_vanleer)
475      call VTb(VThallo)
476     !$OMP END MASTER
[630]477
[1844]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
[630]482
[1844]483      call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm, &
[1549]484             jj_nb_caldyn,0,0,Request_vanleer)
[960]485
[1844]486      call SendRequest(Request_vanleer)
487      !$OMP BARRIER
488      call WaitRequest(Request_vanleer)     
[985]489
[1844]490 !$OMP BARRIER
491 !$OMP MASTER
492      call SetDistrib(jj_nb_caldyn)
493      call VTe(VThallo)
494      call resume_timer(timer_caldyn)
[1549]495 !$OMP END MASTER
496 !$OMP BARRIER 
497        iadvtr=0
498  ENDIF ! if iadvtr.EQ.iapp_tracvl
[630]499
[1549]500END SUBROUTINE advtrac_p
Note: See TracBrowser for help on using the repository browser.