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

Last change on this file since 1750 was 1549, checked in by lguez, 13 years ago

Conversion to free source form, no other change

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