source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/dyn3dpar/advtrac_p.F90 @ 245

Last change on this file since 245 was 222, checked in by ymipsl, 10 years ago

Creating temporary dynamico/lmdz/saturn branche

YM

File size: 17.2 KB
Line 
1! $Id: advtrac_p.F90 1549 2011-07-05 08:41:12Z lguez $
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 instantanee 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 diagnostiques 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!!! ATTENTION !!!! TOUT CE QUI EST ENTRE ICI ET 1234 EST OBSOLETE !!!!!!!
271     GOTO 1234
272!!! ATTENTION !!!!
273
274     !-----------------------------------------------------------
275     !     Appel des sous programmes d'advection
276     !-----------------------------------------------------------
277     do iq=1,nqtot
278        !        call clock(t_initial)
279        if(iadv(iq) == 0) cycle 
280        !   ----------------------------------------------------------------
281        !   Schema de Van Leer I MUSCL
282        !   ----------------------------------------------------------------
283        if(iadv(iq).eq.10) THEN
284
285           call vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
286
287           !   ----------------------------------------------------------------
288           !   Schema "pseudo amont" + test sur humidite specifique
289           !    pour la vapeur d'eau. F. Codron
290           !   ----------------------------------------------------------------
291        else if(iadv(iq).eq.14) then
292           !
293           !ym        stop 'advtrac : appel à vlspltqs :schema non parallelise'
294           CALL vlspltqs_p( q(1,1,1), 2., massem, wg , &
295                pbarug,pbarvg,dtvr,p_tmp,pk_tmp,teta_tmp )
296           !   ----------------------------------------------------------------
297           !   Schema de Frederic Hourdin
298           !   ----------------------------------------------------------------
299        else if(iadv(iq).eq.12) then
300           stop 'advtrac : schema non parallelise'
301           !            Pas de temps adaptatif
302           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
303           if (n.GT.1) then
304              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
305                   dtvr,'n=',n
306           endif
307           do indice=1,n
308              call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
309           end do
310        else if(iadv(iq).eq.13) then
311           stop 'advtrac : schema non parallelise'
312           !            Pas de temps adaptatif
313           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
314           if (n.GT.1) then
315              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
316                   dtvr,'n=',n
317           endif
318           do indice=1,n
319              call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
320           end do
321           !   ----------------------------------------------------------------
322           !   Schema de pente SLOPES
323           !   ----------------------------------------------------------------
324        else if (iadv(iq).eq.20) then
325           stop 'advtrac : schema non parallelise'
326
327           call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
328
329           !   ----------------------------------------------------------------
330           !   Schema de Prather
331           !   ----------------------------------------------------------------
332        else if (iadv(iq).eq.30) then
333           stop 'advtrac : schema non parallelise'
334           !            Pas de temps adaptatif
335           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
336           if (n.GT.1) then
337              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
338                   dtvr,'n=',n
339           endif
340           call  prather(q(1,1,iq),wg,massem,pbarug,pbarvg, &
341                n,dtbon)
342           !   ----------------------------------------------------------------
343           !   Schemas PPM Lin et Rood
344           !   ----------------------------------------------------------------
345        else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND. &
346             iadv(iq).LE.18)) then
347
348           stop 'advtrac : schema non parallelise'
349
350           !        Test sur le flux horizontal
351           !        Pas de temps adaptatif
352           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
353           if (n.GT.1) then
354              write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=', &
355                   dtvr,'n=',n
356           endif
357           !        Test sur le flux vertical
358           CFLmaxz=0.
359           do l=2,llm
360              do ij=iip2,ip1jm
361                 aaa=wg(ij,l)*dtvr/massem(ij,l)
362                 CFLmaxz=max(CFLmaxz,aaa)
363                 bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
364                 CFLmaxz=max(CFLmaxz,bbb)
365              enddo
366           enddo
367           if (CFLmaxz.GE.1) then
368              write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
369           endif
370
371           !-----------------------------------------------------------
372           !        Ss-prg interface LMDZ.4->PPM3d
373           !-----------------------------------------------------------
374
375           call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem, &
376                apppm,bpppm,massebx,masseby,pbarug,pbarvg, &
377                unatppm,vnatppm,psppm)
378
379           do indice=1,n
380              !----------------------------------------------------------------
381              !                         VL (version PPM) horiz. et PPM vert.
382              !----------------------------------------------------------------
383              if (iadv(iq).eq.11) then
384                 !                  Ss-prg PPM3d de Lin
385                 call ppm3d(1,qppm(1,1,iq), &
386                      psppm,psppm, &
387                      unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1, &
388                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
389                      fill,dum,220.)
390
391                 !-------------------------------------------------------------
392                 !                           Monotonic PPM
393                 !-------------------------------------------------------------
394              else if (iadv(iq).eq.16) then
395                 !                  Ss-prg PPM3d de Lin
396                 call ppm3d(1,qppm(1,1,iq), &
397                      psppm,psppm, &
398                      unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1, &
399                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
400                      fill,dum,220.)
401                 !-------------------------------------------------------------
402
403                 !-------------------------------------------------------------
404                 !                           Semi Monotonic PPM
405                 !-------------------------------------------------------------
406              else if (iadv(iq).eq.17) then
407                 !                  Ss-prg PPM3d de Lin
408                 call ppm3d(1,qppm(1,1,iq), &
409                      psppm,psppm, &
410                      unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1, &
411                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
412                      fill,dum,220.)
413                 !-------------------------------------------------------------
414
415                 !-------------------------------------------------------------
416                 !                         Positive Definite PPM
417                 !-------------------------------------------------------------
418              else if (iadv(iq).eq.18) then
419                 !                  Ss-prg PPM3d de Lin
420                 call ppm3d(1,qppm(1,1,iq), &
421                      psppm,psppm, &
422                      unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1, &
423                      iim,jjp1,2,llm,apppm,bpppm,0.01,6400000, &
424                      fill,dum,220.)
425                 !-------------------------------------------------------------
426              endif
427           enddo
428           !-----------------------------------------------------------------
429           !               Ss-prg interface PPM3d-LMDZ.4
430           !-----------------------------------------------------------------
431           call interpost(q(1,1,iq),qppm(1,1,iq))
432        endif
433        !----------------------------------------------------------------------
434
435        !-----------------------------------------------------------------
436        ! On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
437        ! et Nord j=1
438        !-----------------------------------------------------------------
439
440        !                  call traceurpole(q(1,1,iq),massem)
441
442        ! calcul du temps cpu pour un schema donne
443
444        !                  call clock(t_final)
445        !ym                  tps_cpu=t_final-t_initial
446        !ym                  cpuadv(iq)=cpuadv(iq)+tps_cpu
447
448     end DO
449
450!!! ATTENTION !!!!
4511234 CONTINUE
452!!! ATTENTION !!!! LE CODE REPREND ICI !!!!!!!!
453
454     !$OMP BARRIER
455
456     if (planet_type=="earth") then
457
458        ijb=ij_begin
459        ije=ij_end
460
461        !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
462        DO l = 1, llm
463           DO ij = ijb, ije
464              finmasse(ij,l) =  p(ij,l) - p(ij,l+1) 
465           ENDDO
466        ENDDO
467        !$OMP END DO
468
469        CALL qminimum_p( q, 2, finmasse )
470
471     endif ! of if (planet_type=="earth")
472
473        !------------------------------------------------------------------
474        !   on reinitialise a zero les flux de masse cumules
475        !---------------------------------------------------
476        !          iadvtr=0
477
478        !$OMP MASTER
479        call VTe(VTadvection)
480        call stop_timer(timer_vanleer)
481        call VTb(VThallo)
482        !$OMP END MASTER
483
484
485        do j=1,nqtot
486           call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm, &
487                jj_nb_caldyn,0,0,Request_vanleer)
488        enddo
489
490        call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm, &
491             jj_nb_caldyn,0,0,Request_vanleer)
492
493        call SendRequest(Request_vanleer)
494        !$OMP BARRIER
495        call WaitRequest(Request_vanleer)     
496
497        !$OMP BARRIER
498        !$OMP MASTER
499        call SetDistrib(jj_nb_caldyn)
500        call VTe(VThallo)
501        call resume_timer(timer_caldyn)
502 !$OMP END MASTER
503 !$OMP BARRIER 
504        iadvtr=0
505  ENDIF ! if iadvtr.EQ.iapp_tracvl
506
507END SUBROUTINE advtrac_p
Note: See TracBrowser for help on using the repository browser.