source: trunk/LMDZ.COMMON/libf/dyn3dpar/advtrac_p.F90 @ 2236

Last change on this file since 2236 was 2126, checked in by emillour, 6 years ago

Common dynamics:
Some work to enforce total tracer mass conservation in the dynamics.
Still to be further studied and validated.
For now these changes are triggered by setting a "force_conserv_tracer"
flag to ".true." in run.def (default is ".false." to not change anything
with respect to previous versions).
When force_conserv_tracer=.true. then:

  1. Rescale tracer mass in caladvtrac after tracer advection computations
  2. Recompute q ratios once atmospheric mass has been updated in integrd

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