source: LMDZ4/trunk/libf/dyn3dpar/vlspltgen_p.F @ 1119

Last change on this file since 1119 was 985, checked in by Laurent Fairhead, 16 years ago

Mise a jour de dyn3dpar par rapport a dyn3d, inclusion OpenMP et filtre FFT YM
LF
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 10.4 KB
RevLine 
[763]1!
2! $Header$
3!
4       SUBROUTINE vlspltgen_p( q,iadv,pente_max,masse,w,pbaru,pbarv,pdt,
5     ,                                  p,pk,teta                 )
6c
7c     Auteurs:   P.Le Van, F.Hourdin, F.Forget, F.Codron
8c
9c    ********************************************************************
10c          Shema  d'advection " pseudo amont " .
11c      + test sur humidite specifique: Q advecte< Qsat aval
12c                   (F. Codron, 10/99)
13c    ********************************************************************
14c     q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
15c
16c     pente_max facteur de limitation des pentes: 2 en general
17c                                                0 pour un schema amont
18c     pbaru,pbarv,w flux de masse en u ,v ,w
19c     pdt pas de temps
20c
21c     teta temperature potentielle, p pression aux interfaces,
22c     pk exner au milieu des couches necessaire pour calculer Qsat
23c   --------------------------------------------------------------------
24      USE parallel
25      USE mod_hallo
26      USE Write_Field_p
27      USE VAMPIR
28      IMPLICIT NONE
29
30c
31#include "dimensions.h"
32#include "paramet.h"
33#include "logic.h"
34#include "comvert.h"
35#include "comconst.h"
36
37c
38c   Arguments:
39c   ----------
40      INTEGER iadv(nqmx)
41      REAL masse(ip1jmp1,llm),pente_max
42      REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
43      REAL q(ip1jmp1,llm,nqmx)
44      REAL w(ip1jmp1,llm),pdt
45      REAL p(ip1jmp1,llmp1),teta(ip1jmp1,llm),pk(ip1jmp1,llm)
46c
47c      Local
48c   ---------
49c
50      INTEGER ij,l
51c
52      REAL,SAVE :: qsat(ip1jmp1,llm)
53      REAL,SAVE :: zm(ip1jmp1,llm,nqmx)
54      REAL,SAVE :: mu(ip1jmp1,llm)
55      REAL,SAVE :: mv(ip1jm,llm)
56      REAL,SAVE :: mw(ip1jmp1,llm+1)
57      REAL,SAVE :: zq(ip1jmp1,llm,nqmx)
58      REAL zzpbar, zzw
59
60      REAL qmin,qmax
61      DATA qmin,qmax/0.,1.e33/
62
63c--pour rapport de melange saturant--
64
65      REAL rtt,retv,r2es,r3les,r3ies,r4les,r4ies,play
66      REAL ptarg,pdelarg,foeew,zdelta
67      REAL tempe(ip1jmp1)
68      INTEGER ijb,ije,iq
69      type(request) :: MyRequest1
70      type(request) :: MyRequest2
71
72c    fonction psat(T)
73
74       FOEEW ( PTARG,PDELARG ) = EXP (
75     *          (R3LES*(1.-PDELARG)+R3IES*PDELARG) * (PTARG-RTT)
76     * / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG)) )
77
78        r2es  = 380.11733
79        r3les = 17.269
80        r3ies = 21.875
81        r4les = 35.86
82        r4ies = 7.66
83        retv = 0.6077667
84        rtt  = 273.16
85
86c-- Calcul de Qsat en chaque point
87c-- approximation: au milieu des couches play(l)=(p(l)+p(l+1))/2
88c   pour eviter une exponentielle.
[985]89
[763]90      call SetTag(MyRequest1,100)
91      call SetTag(MyRequest2,101)
[985]92
[763]93       
94        ijb=ij_begin-iip1
95        ije=ij_end+iip1
96        if (pole_nord) ijb=ij_begin
97        if (pole_sud) ije=ij_end
98       
99c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
100        DO l = 1, llm
101         DO ij = ijb, ije
102          tempe(ij) = teta(ij,l) * pk(ij,l) /cpp
103         ENDDO
104         DO ij = ijb, ije
105          zdelta = MAX( 0., SIGN(1., rtt - tempe(ij)) )
106          play   = 0.5*(p(ij,l)+p(ij,l+1))
107          qsat(ij,l) = MIN(0.5, r2es* FOEEW(tempe(ij),zdelta) / play )
108          qsat(ij,l) = qsat(ij,l) / ( 1. - retv * qsat(ij,l) )
109         ENDDO
110        ENDDO
111c$OMP END DO NOWAIT
112c      PRINT*,'Debut vlsplt version debug sans vlyqs'
113
114        zzpbar = 0.5 * pdt
115        zzw    = pdt
116
117      ijb=ij_begin
118      ije=ij_end
119      if (pole_nord) ijb=ijb+iip1
120      if (pole_sud)  ije=ije-iip1
121
122c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
123      DO l=1,llm
124        DO ij = ijb,ije
125            mu(ij,l)=pbaru(ij,l) * zzpbar
126         ENDDO
127      ENDDO
128c$OMP END DO NOWAIT
129     
130      ijb=ij_begin-iip1
131      ije=ij_end
132      if (pole_nord) ijb=ij_begin
133      if (pole_sud)  ije=ij_end-iip1
134
135c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
136      DO l=1,llm
137         DO ij=ijb,ije
138            mv(ij,l)=pbarv(ij,l) * zzpbar
139         ENDDO
140      ENDDO
141c$OMP END DO NOWAIT
142
143      ijb=ij_begin
144      ije=ij_end
145
146c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
147      DO l=1,llm
148         DO ij=ijb,ije
149            mw(ij,l)=w(ij,l) * zzw
150         ENDDO
151      ENDDO
152c$OMP END DO NOWAIT
153
154c$OMP MASTER
155      DO ij=ijb,ije
156         mw(ij,llm+1)=0.
157      ENDDO
158c$OMP END MASTER
159
160c      CALL SCOPY(ijp1llm,q,1,zq,1)
161c      CALL SCOPY(ijp1llm,masse,1,zm,1)
162
163       ijb=ij_begin
164       ije=ij_end
165
166      DO iq=1,nqmx
167c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
168        DO l=1,llm
169          zq(ijb:ije,l,iq)=q(ijb:ije,l,iq)
170          zm(ijb:ije,l,iq)=masse(ijb:ije,l)
171        ENDDO
172c$OMP END DO NOWAIT
173      ENDDO
174
175
176c$OMP BARRIER           
177      DO iq=1,nqmx
178
179        if(iadv(iq) == 0) then
180       
181          cycle
182       
183        else if (iadv(iq)==10) then
[985]184
185#ifdef _ADV_HALO       
[763]186          call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
[985]187     &               ij_begin,ij_begin+2*iip1-1)
188          call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
189     &               ij_end-2*iip1+1,ij_end)
190#else
191          call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
[854]192     &               ij_begin,ij_end)
[985]193#endif
[763]194
195c$OMP MASTER
196          call VTb(VTHallo)
[985]197c$OMP END MASTER
[763]198          call Register_Hallo(zq(1,1,iq),ip1jmp1,llm,2,2,2,2,MyRequest1)
199          call Register_Hallo(zm(1,1,iq),ip1jmp1,llm,1,1,1,1,MyRequest1)
[985]200
201c$OMP MASTER
[763]202          call VTe(VTHallo)
203c$OMP END MASTER
204        else if (iadv(iq)==14) then
[985]205
206#ifdef _ADV_HALO           
[763]207          call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat,
[985]208     &                 ij_begin,ij_begin+2*iip1-1)
209          call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat,
210     &                 ij_end-2*iip1+1,ij_end)
211#else
212
213          call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat,
[854]214     &                 ij_begin,ij_end)
[985]215#endif
216
[763]217c$OMP MASTER
218          call VTb(VTHallo)
[985]219c$OMP END MASTER
220
[763]221          call Register_Hallo(zq(1,1,iq),ip1jmp1,llm,2,2,2,2,MyRequest1)
222          call Register_Hallo(zm(1,1,iq),ip1jmp1,llm,1,1,1,1,MyRequest1)
[985]223
224c$OMP MASTER
[763]225          call VTe(VTHallo)
226c$OMP END MASTER
227        else
228       
229          stop 'vlspltgen_p : schema non parallelise'
230     
231        endif
232     
233      enddo
234     
235     
236c$OMP BARRIER     
237c$OMP MASTER     
238      call VTb(VTHallo)
[985]239c$OMP END MASTER
240
[763]241      call SendRequest(MyRequest1)
[985]242
243c$OMP MASTER
[763]244      call VTe(VTHallo)
[985]245c$OMP END MASTER       
246c$OMP BARRIER
247      do iq=1,nqmx
248
249        if(iadv(iq) == 0) then
250       
251          cycle
252       
253        else if (iadv(iq)==10) then
254
255#ifdef _ADV_HALLO
256          call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
257     &               ij_begin+2*iip1,ij_end-2*iip1)
258#endif       
259        else if (iadv(iq)==14) then
260#ifdef _ADV_HALLO
261          call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat,
262     &                 ij_begin+2*iip1,ij_end-2*iip1)
263#endif   
264        else
265       
266          stop 'vlspltgen_p : schema non parallelise'
267     
268        endif
269     
270      enddo
271c$OMP BARRIER     
272c$OMP MASTER
273      call VTb(VTHallo)
[763]274c$OMP END MASTER
[985]275
276!      call WaitRecvRequest(MyRequest1)
277!      call WaitSendRequest(MyRequest1)
[763]278c$OMP BARRIER
[985]279       call WaitRequest(MyRequest1)
280
281
282c$OMP MASTER
283      call VTe(VTHallo)
284c$OMP END MASTER
285c$OMP BARRIER
[763]286 
287      do iq=1,nqmx
288
289        if(iadv(iq) == 0) then
290       
291          cycle
292       
293        else if (iadv(iq)==10) then
294       
295          call vly_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv)
296 
297        else if (iadv(iq)==14) then
298     
299          call vlyqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv,qsat)
300 
301        else
302       
303          stop 'vlspltgen_p : schema non parallelise'
304     
305        endif
306       
307       enddo
308
[985]309
[763]310      do iq=1,nqmx
311
312        if(iadv(iq) == 0) then
313         
314          cycle
315       
316        else if (iadv(iq)==10 .or. iadv(iq)==14 ) then
317
318c$OMP BARRIER       
[985]319#ifdef _ADV_HALLO
[763]320          call vlz_p(zq(1,1,iq),pente_max,zm(1,1,iq),mw,
[985]321     &               ij_begin,ij_begin+2*iip1-1)
322          call vlz_p(zq(1,1,iq),pente_max,zm(1,1,iq),mw,
323     &               ij_end-2*iip1+1,ij_end)
324#else
325          call vlz_p(zq(1,1,iq),pente_max,zm(1,1,iq),mw,
[854]326     &               ij_begin,ij_end)
[985]327#endif
[763]328c$OMP BARRIER
329
330c$OMP MASTER
331          call VTb(VTHallo)
[985]332c$OMP END MASTER
333
[763]334          call Register_Hallo(zq(1,1,iq),ip1jmp1,llm,2,2,2,2,MyRequest2)
335          call Register_Hallo(zm(1,1,iq),ip1jmp1,llm,1,1,1,1,MyRequest2)
[985]336
337c$OMP MASTER
[763]338          call VTe(VTHallo)
339c$OMP END MASTER       
340c$OMP BARRIER
341        else
342       
343          stop 'vlspltgen_p : schema non parallelise'
344     
345        endif
346     
347      enddo
348c$OMP BARRIER     
[985]349
[763]350c$OMP MASTER       
351      call VTb(VTHallo)
[985]352c$OMP END MASTER
353
[763]354      call SendRequest(MyRequest2)
[985]355
356c$OMP MASTER
[763]357      call VTe(VTHallo)
[985]358c$OMP END MASTER       
359
360c$OMP BARRIER
361      do iq=1,nqmx
362
363        if(iadv(iq) == 0) then
364         
365          cycle
366       
367        else if (iadv(iq)==10 .or. iadv(iq)==14 ) then
368c$OMP BARRIER       
369
370#ifdef _ADV_HALLO
371          call vlz_p(zq(1,1,iq),pente_max,zm(1,1,iq),mw,
372     &               ij_begin+2*iip1,ij_end-2*iip1)
373#endif
374
375c$OMP BARRIER       
376        else
377       
378          stop 'vlspltgen_p : schema non parallelise'
379     
380        endif
381     
382      enddo
383
384c$OMP BARRIER
385c$OMP MASTER
386      call VTb(VTHallo)
[763]387c$OMP END MASTER
[985]388
389!      call WaitRecvRequest(MyRequest2)
390!      call WaitSendRequest(MyRequest2)
[763]391c$OMP BARRIER
[985]392       CALL WaitRequest(MyRequest2)
[763]393
[985]394c$OMP MASTER
395      call VTe(VTHallo)
396c$OMP END MASTER
397c$OMP BARRIER
[763]398
[985]399
[763]400      do iq=1,nqmx
401
402        if(iadv(iq) == 0) then
403       
404          cycle
405       
406        else if (iadv(iq)==10) then
407       
408          call vly_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv)
409 
410        else if (iadv(iq)==14) then
411     
412          call vlyqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv,qsat)
413 
414        else
415       
416          stop 'vlspltgen_p : schema non parallelise'
417     
418        endif
419       
420       enddo
421
422      do iq=1,nqmx
423
424        if(iadv(iq) == 0) then
425         
426          cycle
427       
428        else if (iadv(iq)==10) then
429       
430          call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
431     &               ij_begin,ij_end)
432 
433        else if (iadv(iq)==14) then
434     
435          call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat,
436     &                 ij_begin,ij_end)
437 
438        else
439       
440          stop 'vlspltgen_p : schema non parallelise'
441     
442        endif
443       
444       enddo
445
446     
447      ijb=ij_begin
448      ije=ij_end
449c$OMP BARRIER     
450
[985]451
[763]452      DO iq=1,nqmx
453
454c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
455        DO l=1,llm
456           DO ij=ijb,ije
457c             print *,'zq-->',ij,l,iq,zq(ij,l,iq)
458c            print *,'q-->',ij,l,iq,q(ij,l,iq)
459             q(ij,l,iq)=zq(ij,l,iq)
460           ENDDO
461        ENDDO
462c$OMP END DO NOWAIT         
463
464c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
465        DO l=1,llm
466           DO ij=ijb,ije-iip1+1,iip1
467              q(ij+iim,l,iq)=q(ij,l,iq)
468           ENDDO
469        ENDDO
470c$OMP END DO NOWAIT 
471
472      ENDDO
473       
474     
475c$OMP BARRIER
476
477cc$OMP MASTER     
478c      call WaitSendRequest(MyRequest1)
479c      call WaitSendRequest(MyRequest2)
480cc$OMP END MASTER
481cc$OMP BARRIER
482
483      RETURN
484      END
Note: See TracBrowser for help on using the repository browser.