source: LMDZ5/trunk/libf/dyn3dmem/advtrac_loc.F @ 2052

Last change on this file since 2052 was 1987, checked in by Ehouarn Millour, 11 years ago

Add updating pressure, mass and Exner function (ie: all variables which depend on surface pressure) after adding physics tendencies (which include a surface pressure tendency).
Note that this change induces slight changes in GCM results with respect to previous svn version of the code, even if surface pressure tendency is zero (because of recomputation of polar values as an average over polar points on the dynamics grid).
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
File size: 12.7 KB
Line 
1!
2! $Id$
3!
4c
5c
6#define DEBUG_IO
7#undef DEBUG_IO
8      SUBROUTINE advtrac_loc(pbarug,pbarvg ,wg,
9     *                   p,  massem,q,teta,
10     *                   pk   )
11
12c     Auteur :  F. Hourdin
13c
14c     Modif. P. Le Van     (20/12/97)
15c            F. Codron     (10/99)
16c            D. Le Croller (07/2001)
17c            M.A Filiberti (04/2002)
18c
19      USE parallel_lmdz
20      USE Write_Field_loc
21      USE Write_Field
22      USE Bands
23      USE mod_hallo
24      USE Vampir
25      USE times
26      USE infotrac, ONLY: nqtot, iadv
27      USE control_mod, ONLY: iapp_tracvl, day_step, planet_type
28      USE advtrac_mod, ONLY: finmasse
29      IMPLICIT NONE
30c
31#include "dimensions.h"
32#include "paramet.h"
33#include "comconst.h"
34#include "comvert.h"
35#include "comdissip.h"
36#include "comgeom2.h"
37#include "logic.h"
38#include "temps.h"
39#include "ener.h"
40#include "description.h"
41
42c-------------------------------------------------------------------
43c     Arguments
44c-------------------------------------------------------------------
45c     Ajout PPM
46c--------------------------------------------------------
47      REAL massebx(ijb_u:ije_u,llm),masseby(ijb_v:ije_v,llm)
48c--------------------------------------------------------
49      INTEGER iapptrac
50      REAL pbarug(ijb_u:ije_u,llm),pbarvg(ijb_v:ije_v,llm)
51      REAL wg(ijb_u:ije_u,llm)
52      REAL q(ijb_u:ije_u,llm,nqtot),massem(ijb_u:ije_u,llm)
53      REAL p( ijb_u:ije_u,llmp1 ),teta(ijb_u:ije_u,llm)
54      REAL pk(ijb_u:ije_u,llm)
55
56c-------------------------------------------------------------
57c     Variables locales
58c-------------------------------------------------------------
59
60      REAL zdp(ijb_u:ije_u)
61      REAL (kind=kind(1.d0)) :: t_initial, t_final, tps_cpu
62      INTEGER,SAVE :: iadvtr=0
63c$OMP THREADPRIVATE(iadvtr)
64      INTEGER ij,l,iq,iiq
65      REAL zdpmin, zdpmax
66c----------------------------------------------------------
67c     Rajouts pour PPM
68c----------------------------------------------------------
69      INTEGER indice,n
70      REAL dtbon ! Pas de temps adaptatif pour que CFL<1
71      REAL CFLmaxz,aaa,bbb ! CFL maximum
72      REAL psppm(iim,jjb_u:jje_u) ! pression  au sol
73      REAL unatppm(iim,jjb_u:jje_u,llm),vnatppm(iim,jjb_u:jje_u,llm)
74      REAL qppm(iim*jjnb_u,llm,nqtot)
75      REAL fluxwppm(iim,jjb_u:jje_u,llm)
76      REAL apppm(llmp1), bpppm(llmp1)
77      LOGICAL dum,fill
78      DATA fill/.true./
79      DATA dum/.true./
80      integer ijb,ije,ijbu,ijbv,ijeu,ijev,j
81      type(Request),SAVE :: testRequest
82!$OMP THREADPRIVATE(testRequest)
83
84c  test sur l'eventuelle creation de valeurs negatives de la masse
85         ijb=ij_begin
86         ije=ij_end
87         if (pole_nord) ijb=ij_begin+iip1
88         if (pole_sud) ije=ij_end-iip1
89         
90c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
91         DO l=1,llm-1
92            DO ij = ijb+1,ije
93              zdp(ij) =    pbarug(ij-1,l)   - pbarug(ij,l)
94     s                  - pbarvg(ij-iip1,l) + pbarvg(ij,l)
95     s                  +       wg(ij,l+1)  - wg(ij,l)
96            ENDDO
97           
98c            CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
99c ym  ---> pourquoi jjm-1 et non jjm ? a cause du pole ?
100           
101            do ij=ijb,ije-iip1+1,iip1
102              zdp(ij)=zdp(ij+iip1-1)
103            enddo
104           
105            DO ij = ijb,ije
106               zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
107            ENDDO
108
109
110c            CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
111c  ym ---> eventuellement a revoir
112            CALL minmax ( ije-ijb+1, zdp(ijb), zdpmin,zdpmax )
113           
114            IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN
115            PRINT*,'WARNING DP/P l=',l,'  MIN:',zdpmin,
116     s        '   MAX:', zdpmax
117            ENDIF
118
119         ENDDO
120c$OMP END DO NOWAIT
121
122c-------------------------------------------------------------------
123c   Advection proprement dite (Modification Le Croller (07/2001)
124c-------------------------------------------------------------------
125
126c----------------------------------------------------
127c        Calcul des moyennes basées sur la masse
128c----------------------------------------------------
129
130cym      ----> Normalement, inutile pour les schémas classiques
131cym      ----> Revérifier lors de la parallélisation des autres schemas
132   
133cym          call massbar_p(massem,massebx,masseby)         
134
135#ifdef DEBUG_IO   
136          CALL WriteField_u('massem',massem)
137          CALL WriteField_u('wg',wg)
138          CALL WriteField_u('pbarug',pbarug)
139          CALL WriteField_v('pbarvg',pbarvg)
140          CALL WriteField_u('p_tmp',p)
141          CALL WriteField_u('pk_tmp',pk)
142          CALL WriteField_u('teta_tmp',teta)
143          do j=1,nqtot
144            call WriteField_u('q_adv'//trim(int2str(j)),
145     .                q(:,:,j))
146          enddo
147#endif
148
149!         
150!       call Register_Hallo_v(pbarvg,llm,1,1,1,1,TestRequest)
151!
152!       call SendRequest(TestRequest)
153!c$OMP BARRIER
154!       call WaitRequest(TestRequest)
155c$OMP BARRIER
156                 
157          call vlspltgen_loc( q,iadv, 2., massem, wg ,
158     *                        pbarug,pbarvg,dtvr,p,
159     *                        pk,teta )
160
161#ifdef DEBUG_IO     
162          do j=1,nqtot
163            call WriteField_u('q_adv'//trim(int2str(j)),
164     .                q(:,:,j))
165          enddo
166#endif
167         
168          GOTO 1234     
169c-----------------------------------------------------------
170c     Appel des sous programmes d'advection
171c-----------------------------------------------------------
172      do iq=1,nqtot
173c        call clock(t_initial)
174        if(iadv(iq) == 0) cycle
175c   ----------------------------------------------------------------
176c   Schema de Van Leer I MUSCL
177c   ----------------------------------------------------------------
178        if(iadv(iq).eq.10) THEN
179     
180!LF         call vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
181
182c   ----------------------------------------------------------------
183c   Schema "pseudo amont" + test sur humidite specifique
184C    pour la vapeur d'eau. F. Codron
185c   ----------------------------------------------------------------
186        else if(iadv(iq).eq.14) then
187c
188cym           stop 'advtrac : appel à vlspltqs :schema non parallelise'
189!LF           CALL vlspltqs_p( q(1,1,1), 2., massem, wg ,
190!LF     *                 pbarug,pbarvg,dtvr,p,pk,teta )
191c   ----------------------------------------------------------------
192c   Schema de Frederic Hourdin
193c   ----------------------------------------------------------------
194        else if(iadv(iq).eq.12) then
195          stop 'advtrac : schema non parallelise'
196c            Pas de temps adaptatif
197           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
198           if (n.GT.1) then
199           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
200     s             dtvr,'n=',n
201           endif
202           do indice=1,n
203            call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
204           end do
205        else if(iadv(iq).eq.13) then
206          stop 'advtrac : schema non parallelise'
207c            Pas de temps adaptatif
208           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
209           if (n.GT.1) then
210           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
211     s             dtvr,'n=',n
212           endif
213          do indice=1,n
214            call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
215          end do
216c   ----------------------------------------------------------------
217c   Schema de pente SLOPES
218c   ----------------------------------------------------------------
219        else if (iadv(iq).eq.20) then
220          stop 'advtrac : schema non parallelise'
221
222            call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
223
224c   ----------------------------------------------------------------
225c   Schema de Prather
226c   ----------------------------------------------------------------
227        else if (iadv(iq).eq.30) then
228          stop 'advtrac : schema non parallelise'
229c            Pas de temps adaptatif
230           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
231           if (n.GT.1) then
232           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
233     s             dtvr,'n=',n
234           endif
235           call  prather(q(1,1,iq),wg,massem,pbarug,pbarvg,
236     s                     n,dtbon)
237c   ----------------------------------------------------------------
238c   Schemas PPM Lin et Rood
239c   ----------------------------------------------------------------
240       else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND.
241     s                     iadv(iq).LE.18)) then
242
243           stop 'advtrac : schema non parallelise'
244
245c        Test sur le flux horizontal
246c        Pas de temps adaptatif
247         call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
248         if (n.GT.1) then
249           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
250     s             dtvr,'n=',n
251         endif
252c        Test sur le flux vertical
253         CFLmaxz=0.
254         do l=2,llm
255           do ij=iip2,ip1jm
256            aaa=wg(ij,l)*dtvr/massem(ij,l)
257            CFLmaxz=max(CFLmaxz,aaa)
258            bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
259            CFLmaxz=max(CFLmaxz,bbb)
260           enddo
261         enddo
262         if (CFLmaxz.GE.1) then
263            write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
264         endif
265
266c-----------------------------------------------------------
267c        Ss-prg interface LMDZ.4->PPM3d
268c-----------------------------------------------------------
269
270          call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem,
271     s                 apppm,bpppm,massebx,masseby,pbarug,pbarvg,
272     s                 unatppm,vnatppm,psppm)
273
274          do indice=1,n
275c---------------------------------------------------------------------
276c                         VL (version PPM) horiz. et PPM vert.
277c---------------------------------------------------------------------
278                if (iadv(iq).eq.11) then
279c                  Ss-prg PPM3d de Lin
280                  call ppm3d(1,qppm(1,1,iq),
281     s                       psppm,psppm,
282     s                       unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1,
283     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
284     s                       fill,dum,220.)
285
286c----------------------------------------------------------------------
287c                           Monotonic PPM
288c----------------------------------------------------------------------
289               else if (iadv(iq).eq.16) then
290c                  Ss-prg PPM3d de Lin
291                  call ppm3d(1,qppm(1,1,iq),
292     s                       psppm,psppm,
293     s                       unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1,
294     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
295     s                       fill,dum,220.)
296c---------------------------------------------------------------------
297
298c---------------------------------------------------------------------
299c                           Semi Monotonic PPM
300c---------------------------------------------------------------------
301               else if (iadv(iq).eq.17) then
302c                  Ss-prg PPM3d de Lin
303                  call ppm3d(1,qppm(1,1,iq),
304     s                       psppm,psppm,
305     s                       unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1,
306     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
307     s                       fill,dum,220.)
308c---------------------------------------------------------------------
309
310c---------------------------------------------------------------------
311c                         Positive Definite PPM
312c---------------------------------------------------------------------
313                else if (iadv(iq).eq.18) then
314c                  Ss-prg PPM3d de Lin
315                  call ppm3d(1,qppm(1,1,iq),
316     s                       psppm,psppm,
317     s                       unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1,
318     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
319     s                       fill,dum,220.)
320c---------------------------------------------------------------------
321                endif
322            enddo
323c-----------------------------------------------------------------
324c               Ss-prg interface PPM3d-LMDZ.4
325c-----------------------------------------------------------------
326                  call interpost(q(1,1,iq),qppm(1,1,iq))
327            endif
328c----------------------------------------------------------------------
329
330c-----------------------------------------------------------------
331c On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
332c et Nord j=1
333c-----------------------------------------------------------------
334
335c                  call traceurpole(q(1,1,iq),massem)
336
337c calcul du temps cpu pour un schema donne
338
339
340       end DO
341
3421234  CONTINUE
343c$OMP BARRIER
344
345      if (planet_type=="earth") then
346
347      ijb=ij_begin
348      ije=ij_end
349
350c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
351       DO l = 1, llm
352         DO ij = ijb, ije
353           finmasse(ij,l) =  p(ij,l) - p(ij,l+1)
354         ENDDO
355       ENDDO
356c$OMP END DO
357
358       CALL qminimum_loc( q, 2, finmasse )
359
360      endif ! of if (planet_type=="earth")
361
362       RETURN
363       END
364
Note: See TracBrowser for help on using the repository browser.