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

Last change on this file since 5308 was 2622, checked in by Ehouarn Millour, 8 years ago

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