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

Last change on this file since 1660 was 1632, checked in by Laurent Fairhead, 12 years ago

Import initial du répertoire dyn3dmem

Attention! ceci n'est qu'une version préliminaire du code "basse mémoire":
le code contenu dans ce répertoire est basé sur la r1320 et a donc besoin
d'être mis à jour par rapport à la dynamique parallèle d'aujourd'hui.
Ce code est toutefois mis à disposition pour circonvenir à des problèmes
de mémoire que certaines configurations du modèle pourraient rencontrer.
Dans l'état, il compile et tourne sur vargas et au CCRT

Initial import of dyn3dmem

Warning! this is just a preliminary version of the memory light code:
it is based on r1320 of the code and thus needs to be updated before
it can replace the present dyn3dpar code. It is nevertheless put at your
disposal to circumvent some memory problems some LMDZ configurations may
encounter. In its present state, it will compile and run on vargas and CCRT

File size: 12.6 KB
RevLine 
[1632]1!
2! $Id: advtrac_p.F 1299 2010-01-20 14:27:21Z fairhead $
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
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
27      USE control_mod
28      USE advtrac_mod
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) :: testRequest
82
83c  test sur l'eventuelle creation de valeurs negatives de la masse
84         ijb=ij_begin
85         ije=ij_end
86         if (pole_nord) ijb=ij_begin+iip1
87         if (pole_sud) ije=ij_end-iip1
88         
89c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
90         DO l=1,llm-1
91            DO ij = ijb+1,ije
92              zdp(ij) =    pbarug(ij-1,l)   - pbarug(ij,l)
93     s                  - pbarvg(ij-iip1,l) + pbarvg(ij,l)
94     s                  +       wg(ij,l+1)  - wg(ij,l)
95            ENDDO
96           
97c            CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
98c ym  ---> pourquoi jjm-1 et non jjm ? a cause du pole ?
99           
100            do ij=ijb,ije-iip1+1,iip1
101              zdp(ij)=zdp(ij+iip1-1)
102            enddo
103           
104            DO ij = ijb,ije
105               zdp(ij)= zdp(ij)*dtvr/ massem(ij,l)
106            ENDDO
107
108
109c            CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
110c  ym ---> eventuellement a revoir
111            CALL minmax ( ije-ijb+1, zdp(ijb), zdpmin,zdpmax )
112           
113            IF(MAX(ABS(zdpmin),ABS(zdpmax)).GT.0.5) THEN
114            PRINT*,'WARNING DP/P l=',l,'  MIN:',zdpmin,
115     s        '   MAX:', zdpmax
116            ENDIF
117
118         ENDDO
119c$OMP END DO NOWAIT
120
121c-------------------------------------------------------------------
122c   Advection proprement dite (Modification Le Croller (07/2001)
123c-------------------------------------------------------------------
124
125c----------------------------------------------------
126c        Calcul des moyennes basées sur la masse
127c----------------------------------------------------
128
129cym      ----> Normalement, inutile pour les schémas classiques
130cym      ----> Revérifier lors de la parallélisation des autres schemas
131   
132cym          call massbar_p(massem,massebx,masseby)         
133
134#ifdef DEBUG_IO   
135          CALL WriteField_u('massem',massem)
136          CALL WriteField_u('wg',wg)
137          CALL WriteField_u('pbarug',pbarug)
138          CALL WriteField_v('pbarvg',pbarvg)
139          CALL WriteField_u('p_tmp',p)
140          CALL WriteField_u('pk_tmp',pk)
141          CALL WriteField_u('teta_tmp',teta)
142          do j=1,nqtot
143            call WriteField_u('q_adv'//trim(int2str(j)),
144     .                q(:,:,j))
145          enddo
146#endif
147
148!         
149!       call Register_Hallo_v(pbarvg,llm,1,1,1,1,TestRequest)
150!
151!       call SendRequest(TestRequest)
152!c$OMP BARRIER
153!       call WaitRequest(TestRequest)
154c$OMP BARRIER
155                 
156          call vlspltgen_loc( q,iadv, 2., massem, wg ,
157     *                        pbarug,pbarvg,dtvr,p,
158     *                        pk,teta )
159
160#ifdef DEBUG_IO     
161          do j=1,nqtot
162            call WriteField_u('q_adv'//trim(int2str(j)),
163     .                q(:,:,j))
164          enddo
165#endif
166         
167          GOTO 1234     
168c-----------------------------------------------------------
169c     Appel des sous programmes d'advection
170c-----------------------------------------------------------
171      do iq=1,nqtot
172c        call clock(t_initial)
173        if(iadv(iq) == 0) cycle
174c   ----------------------------------------------------------------
175c   Schema de Van Leer I MUSCL
176c   ----------------------------------------------------------------
177        if(iadv(iq).eq.10) THEN
178     
179!LF         call vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
180
181c   ----------------------------------------------------------------
182c   Schema "pseudo amont" + test sur humidite specifique
183C    pour la vapeur d'eau. F. Codron
184c   ----------------------------------------------------------------
185        else if(iadv(iq).eq.14) then
186c
187cym           stop 'advtrac : appel à vlspltqs :schema non parallelise'
188!LF           CALL vlspltqs_p( q(1,1,1), 2., massem, wg ,
189!LF     *                 pbarug,pbarvg,dtvr,p,pk,teta )
190c   ----------------------------------------------------------------
191c   Schema de Frederic Hourdin
192c   ----------------------------------------------------------------
193        else if(iadv(iq).eq.12) then
194          stop 'advtrac : schema non parallelise'
195c            Pas de temps adaptatif
196           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
197           if (n.GT.1) then
198           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
199     s             dtvr,'n=',n
200           endif
201           do indice=1,n
202            call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,1)
203           end do
204        else if(iadv(iq).eq.13) then
205          stop 'advtrac : schema non parallelise'
206c            Pas de temps adaptatif
207           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
208           if (n.GT.1) then
209           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
210     s             dtvr,'n=',n
211           endif
212          do indice=1,n
213            call advn(q(1,1,iq),massem,wg,pbarug,pbarvg,dtbon,2)
214          end do
215c   ----------------------------------------------------------------
216c   Schema de pente SLOPES
217c   ----------------------------------------------------------------
218        else if (iadv(iq).eq.20) then
219          stop 'advtrac : schema non parallelise'
220
221            call pentes_ini (q(1,1,iq),wg,massem,pbarug,pbarvg,0)
222
223c   ----------------------------------------------------------------
224c   Schema de Prather
225c   ----------------------------------------------------------------
226        else if (iadv(iq).eq.30) then
227          stop 'advtrac : schema non parallelise'
228c            Pas de temps adaptatif
229           call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
230           if (n.GT.1) then
231           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
232     s             dtvr,'n=',n
233           endif
234           call  prather(q(1,1,iq),wg,massem,pbarug,pbarvg,
235     s                     n,dtbon)
236c   ----------------------------------------------------------------
237c   Schemas PPM Lin et Rood
238c   ----------------------------------------------------------------
239       else if (iadv(iq).eq.11.OR.(iadv(iq).GE.16.AND.
240     s                     iadv(iq).LE.18)) then
241
242           stop 'advtrac : schema non parallelise'
243
244c        Test sur le flux horizontal
245c        Pas de temps adaptatif
246         call adaptdt(iadv(iq),dtbon,n,pbarug,massem)
247         if (n.GT.1) then
248           write(*,*) 'WARNING horizontal dt=',dtbon,'dtvr=',
249     s             dtvr,'n=',n
250         endif
251c        Test sur le flux vertical
252         CFLmaxz=0.
253         do l=2,llm
254           do ij=iip2,ip1jm
255            aaa=wg(ij,l)*dtvr/massem(ij,l)
256            CFLmaxz=max(CFLmaxz,aaa)
257            bbb=-wg(ij,l)*dtvr/massem(ij,l-1)
258            CFLmaxz=max(CFLmaxz,bbb)
259           enddo
260         enddo
261         if (CFLmaxz.GE.1) then
262            write(*,*) 'WARNING vertical','CFLmaxz=', CFLmaxz
263         endif
264
265c-----------------------------------------------------------
266c        Ss-prg interface LMDZ.4->PPM3d
267c-----------------------------------------------------------
268
269          call interpre(q(1,1,iq),qppm(1,1,iq),wg,fluxwppm,massem,
270     s                 apppm,bpppm,massebx,masseby,pbarug,pbarvg,
271     s                 unatppm,vnatppm,psppm)
272
273          do indice=1,n
274c---------------------------------------------------------------------
275c                         VL (version PPM) horiz. et PPM vert.
276c---------------------------------------------------------------------
277                if (iadv(iq).eq.11) then
278c                  Ss-prg PPM3d de Lin
279                  call ppm3d(1,qppm(1,1,iq),
280     s                       psppm,psppm,
281     s                       unatppm,vnatppm,fluxwppm,dtbon,2,2,2,1,
282     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
283     s                       fill,dum,220.)
284
285c----------------------------------------------------------------------
286c                           Monotonic PPM
287c----------------------------------------------------------------------
288               else if (iadv(iq).eq.16) then
289c                  Ss-prg PPM3d de Lin
290                  call ppm3d(1,qppm(1,1,iq),
291     s                       psppm,psppm,
292     s                       unatppm,vnatppm,fluxwppm,dtbon,3,3,3,1,
293     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
294     s                       fill,dum,220.)
295c---------------------------------------------------------------------
296
297c---------------------------------------------------------------------
298c                           Semi Monotonic PPM
299c---------------------------------------------------------------------
300               else if (iadv(iq).eq.17) then
301c                  Ss-prg PPM3d de Lin
302                  call ppm3d(1,qppm(1,1,iq),
303     s                       psppm,psppm,
304     s                       unatppm,vnatppm,fluxwppm,dtbon,4,4,4,1,
305     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
306     s                       fill,dum,220.)
307c---------------------------------------------------------------------
308
309c---------------------------------------------------------------------
310c                         Positive Definite PPM
311c---------------------------------------------------------------------
312                else if (iadv(iq).eq.18) then
313c                  Ss-prg PPM3d de Lin
314                  call ppm3d(1,qppm(1,1,iq),
315     s                       psppm,psppm,
316     s                       unatppm,vnatppm,fluxwppm,dtbon,5,5,5,1,
317     s                       iim,jjp1,2,llm,apppm,bpppm,0.01,6400000,
318     s                       fill,dum,220.)
319c---------------------------------------------------------------------
320                endif
321            enddo
322c-----------------------------------------------------------------
323c               Ss-prg interface PPM3d-LMDZ.4
324c-----------------------------------------------------------------
325                  call interpost(q(1,1,iq),qppm(1,1,iq))
326            endif
327c----------------------------------------------------------------------
328
329c-----------------------------------------------------------------
330c On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
331c et Nord j=1
332c-----------------------------------------------------------------
333
334c                  call traceurpole(q(1,1,iq),massem)
335
336c calcul du temps cpu pour un schema donne
337
338
339       end DO
340
3411234  CONTINUE
342c$OMP BARRIER
343
344      ijb=ij_begin
345      ije=ij_end
346
347c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
348       DO l = 1, llm
349         DO ij = ijb, ije
350           finmasse(ij,l) =  p(ij,l) - p(ij,l+1)
351         ENDDO
352       ENDDO
353c$OMP END DO
354
355       CALL qminimum_loc( q, 2, finmasse )
356
357
358       RETURN
359       END
360
Note: See TracBrowser for help on using the repository browser.