My Project
 All Classes Files Functions Variables Macros
advtrac_loc.F
Go to the documentation of this file.
1 !
2 ! $Id$
3 !
4 c
5 c
6 #define DEBUG_IO
7 #undef DEBUG_IO
8  SUBROUTINE advtrac_loc(pbarug,pbarvg ,wg,
9  * p, massem,q,teta,
10  * pk )
11 
12 c Auteur : F. Hourdin
13 c
14 c Modif. P. Le Van (20/12/97)
15 c F. Codron (10/99)
16 c D. Le Croller (07/2001)
17 c M.A Filiberti (04/2002)
18 c
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
30 c
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 
42 c-------------------------------------------------------------------
43 c Arguments
44 c-------------------------------------------------------------------
45 c Ajout PPM
46 c--------------------------------------------------------
47  REAL massebx(ijb_u:ije_u,llm),masseby(ijb_v:ije_v,llm)
48 c--------------------------------------------------------
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 
56 c-------------------------------------------------------------
57 c Variables locales
58 c-------------------------------------------------------------
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
63 c$OMP THREADPRIVATE(iadvtr)
64  INTEGER ij,l,iq,iiq
65  REAL zdpmin, zdpmax
66 c----------------------------------------------------------
67 c Rajouts pour PPM
68 c----------------------------------------------------------
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 
83 c 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 
89 c$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 
97 c CALL SCOPY( jjm -1 ,zdp(iip1+iip1),iip1,zdp(iip2),iip1 )
98 c 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 
109 c CALL minmax ( ip1jm-iip1, zdp(iip2), zdpmin,zdpmax )
110 c 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
119 c$OMP END DO NOWAIT
120 
121 c-------------------------------------------------------------------
122 c Advection proprement dite (Modification Le Croller (07/2001)
123 c-------------------------------------------------------------------
124 
125 c----------------------------------------------------
126 c Calcul des moyennes basées sur la masse
127 c----------------------------------------------------
128 
129 cym ----> Normalement, inutile pour les schémas classiques
130 cym ----> Revérifier lors de la parallélisation des autres schemas
131 
132 cym 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)
154 c$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
168 c-----------------------------------------------------------
169 c Appel des sous programmes d'advection
170 c-----------------------------------------------------------
171  do iq=1,nqtot
172 c call clock(t_initial)
173  if(iadv(iq) == 0) cycle
174 c ----------------------------------------------------------------
175 c Schema de Van Leer I MUSCL
176 c ----------------------------------------------------------------
177  if(iadv(iq).eq.10) THEN
178 
179 !LF call vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
180 
181 c ----------------------------------------------------------------
182 c Schema "pseudo amont" + test sur humidite specifique
183 C pour la vapeur d'eau. F. Codron
184 c ----------------------------------------------------------------
185  else if(iadv(iq).eq.14) then
186 c
187 cym 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 )
190 c ----------------------------------------------------------------
191 c Schema de Frederic Hourdin
192 c ----------------------------------------------------------------
193  else if(iadv(iq).eq.12) then
194  stop 'advtrac : schema non parallelise'
195 c 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'
206 c 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
215 c ----------------------------------------------------------------
216 c Schema de pente SLOPES
217 c ----------------------------------------------------------------
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 
223 c ----------------------------------------------------------------
224 c Schema de Prather
225 c ----------------------------------------------------------------
226  else if (iadv(iq).eq.30) then
227  stop 'advtrac : schema non parallelise'
228 c 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)
236 c ----------------------------------------------------------------
237 c Schemas PPM Lin et Rood
238 c ----------------------------------------------------------------
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 
244 c Test sur le flux horizontal
245 c 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
251 c 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 
265 c-----------------------------------------------------------
266 c Ss-prg interface LMDZ.4->PPM3d
267 c-----------------------------------------------------------
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
274 c---------------------------------------------------------------------
275 c VL (version PPM) horiz. et PPM vert.
276 c---------------------------------------------------------------------
277  if (iadv(iq).eq.11) then
278 c 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 
285 c----------------------------------------------------------------------
286 c Monotonic PPM
287 c----------------------------------------------------------------------
288  else if (iadv(iq).eq.16) then
289 c 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.)
295 c---------------------------------------------------------------------
296 
297 c---------------------------------------------------------------------
298 c Semi Monotonic PPM
299 c---------------------------------------------------------------------
300  else if (iadv(iq).eq.17) then
301 c 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.)
307 c---------------------------------------------------------------------
308 
309 c---------------------------------------------------------------------
310 c Positive Definite PPM
311 c---------------------------------------------------------------------
312  else if (iadv(iq).eq.18) then
313 c 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.)
319 c---------------------------------------------------------------------
320  endif
321  enddo
322 c-----------------------------------------------------------------
323 c Ss-prg interface PPM3d-LMDZ.4
324 c-----------------------------------------------------------------
325  call interpost(q(1,1,iq),qppm(1,1,iq))
326  endif
327 c----------------------------------------------------------------------
328 
329 c-----------------------------------------------------------------
330 c On impose une seule valeur du traceur au pôle Sud j=jjm+1=jjp1
331 c et Nord j=1
332 c-----------------------------------------------------------------
333 
334 c call traceurpole(q(1,1,iq),massem)
335 
336 c calcul du temps cpu pour un schema donne
337 
338 
339  end DO
340 
341 1234 CONTINUE
342 c$OMP BARRIER
343 
344  if (planet_type=="earth") then
345 
346  ijb=ij_begin
347  ije=ij_end
348 
349 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
350  DO l = 1, llm
351  DO ij = ijb, ije
352  finmasse(ij,l) = p(ij,l) - p(ij,l+1)
353  ENDDO
354  ENDDO
355 c$OMP END DO
356 
357  CALL qminimum_loc( q, 2, finmasse )
358 
359  endif ! of if (planet_type=="earth")
360 
361  RETURN
362  END
363