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

Last change on this file since 825 was 774, checked in by Laurent Fairhead, 18 years ago

Suite du merge entre la version et la HEAD: quelques modifications de
Yann sur le
LF
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 9.6 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.
89c$OMP MASTER
90      call SetTag(MyRequest1,100)
91      call SetTag(MyRequest2,101)
92c$OMP END MASTER
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
184       
185          call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
186     &               ij_begin,ij_begin+2*iip1-1)
187          call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
188     &               ij_end-2*iip1+1,ij_end)
189
190c$OMP MASTER
191          call VTb(VTHallo)
192          call Register_Hallo(zq(1,1,iq),ip1jmp1,llm,2,2,2,2,MyRequest1)
193          call Register_Hallo(zm(1,1,iq),ip1jmp1,llm,1,1,1,1,MyRequest1)
194          call VTe(VTHallo)
195c$OMP END MASTER
196        else if (iadv(iq)==14) then
197     
198          call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat,
199     &                 ij_begin,ij_begin+2*iip1-1)
200          call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat,
201     &                 ij_end-2*iip1+1,ij_end)
202c$OMP MASTER
203          call VTb(VTHallo)
204          call Register_Hallo(zq(1,1,iq),ip1jmp1,llm,2,2,2,2,MyRequest1)
205          call Register_Hallo(zm(1,1,iq),ip1jmp1,llm,1,1,1,1,MyRequest1)
206          call VTe(VTHallo)
207c$OMP END MASTER
208        else
209       
210          stop 'vlspltgen_p : schema non parallelise'
211     
212        endif
213     
214      enddo
215     
216     
217c$OMP BARRIER     
218c$OMP MASTER     
219      call VTb(VTHallo)
220      call SendRequest(MyRequest1)
221      call VTe(VTHallo)
222c$OMP END MASTER       
223c$OMP BARRIER
224      do iq=1,nqmx
225
226        if(iadv(iq) == 0) then
227       
228          cycle
229       
230        else if (iadv(iq)==10) then
231
232          call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
233     &               ij_begin+2*iip1,ij_end-2*iip1)
234       
235        else if (iadv(iq)==14) then
236
237          call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat,
238     &                 ij_begin+2*iip1,ij_end-2*iip1)
239   
240        else
241       
242          stop 'vlspltgen_p : schema non parallelise'
243     
244        endif
245     
246      enddo
247c$OMP BARRIER     
248c$OMP MASTER
249      call VTb(VTHallo)
250      call WaitRecvRequest(MyRequest1)
251      call WaitSendRequest(MyRequest1)
252      call VTe(VTHallo)
253c$OMP END MASTER
254c$OMP BARRIER
255 
256      do iq=1,nqmx
257
258        if(iadv(iq) == 0) then
259       
260          cycle
261       
262        else if (iadv(iq)==10) then
263       
264          call vly_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv)
265 
266        else if (iadv(iq)==14) then
267     
268          call vlyqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv,qsat)
269 
270        else
271       
272          stop 'vlspltgen_p : schema non parallelise'
273     
274        endif
275       
276       enddo
277
278      do iq=1,nqmx
279
280        if(iadv(iq) == 0) then
281         
282          cycle
283       
284        else if (iadv(iq)==10 .or. iadv(iq)==14 ) then
285
286c$OMP BARRIER       
287          call vlz_p(zq(1,1,iq),pente_max,zm(1,1,iq),mw,
288     &               ij_begin,ij_begin+2*iip1-1)
289          call vlz_p(zq(1,1,iq),pente_max,zm(1,1,iq),mw,
290     &               ij_end-2*iip1+1,ij_end)
291c$OMP BARRIER
292
293c$OMP MASTER
294          call VTb(VTHallo)
295          call Register_Hallo(zq(1,1,iq),ip1jmp1,llm,2,2,2,2,MyRequest2)
296          call Register_Hallo(zm(1,1,iq),ip1jmp1,llm,1,1,1,1,MyRequest2)
297          call VTe(VTHallo)
298c$OMP END MASTER       
299c$OMP BARRIER
300        else
301       
302          stop 'vlspltgen_p : schema non parallelise'
303     
304        endif
305     
306      enddo
307c$OMP BARRIER     
308c$OMP MASTER       
309      call VTb(VTHallo)
310      call SendRequest(MyRequest2)
311      call VTe(VTHallo)
312c$OMP END MASTER       
313c$OMP BARRIER
314      do iq=1,nqmx
315
316        if(iadv(iq) == 0) then
317         
318          cycle
319       
320        else if (iadv(iq)==10 .or. iadv(iq)==14 ) then
321c$OMP BARRIER       
322          call vlz_p(zq(1,1,iq),pente_max,zm(1,1,iq),mw,
323     &               ij_begin+2*iip1,ij_end-2*iip1)
324c$OMP BARRIER       
325        else
326       
327          stop 'vlspltgen_p : schema non parallelise'
328     
329        endif
330     
331      enddo
332
333c$OMP BARRIER
334c$OMP MASTER
335      call VTb(VTHallo)
336      call WaitRecvRequest(MyRequest2)
337      call WaitSendRequest(MyRequest2)
338      call VTe(VTHallo)
339c$OMP END MASTER
340c$OMP BARRIER
341
342
343      do iq=1,nqmx
344
345        if(iadv(iq) == 0) then
346       
347          cycle
348       
349        else if (iadv(iq)==10) then
350       
351          call vly_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv)
352 
353        else if (iadv(iq)==14) then
354     
355          call vlyqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv,qsat)
356 
357        else
358       
359          stop 'vlspltgen_p : schema non parallelise'
360     
361        endif
362       
363       enddo
364
365      do iq=1,nqmx
366
367        if(iadv(iq) == 0) then
368         
369          cycle
370       
371        else if (iadv(iq)==10) then
372       
373          call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
374     &               ij_begin,ij_end)
375 
376        else if (iadv(iq)==14) then
377     
378          call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat,
379     &                 ij_begin,ij_end)
380 
381        else
382       
383          stop 'vlspltgen_p : schema non parallelise'
384     
385        endif
386       
387       enddo
388
389     
390      ijb=ij_begin
391      ije=ij_end
392c$OMP BARRIER     
393
394      DO iq=1,nqmx
395
396c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
397        DO l=1,llm
398           DO ij=ijb,ije
399c             print *,'zq-->',ij,l,iq,zq(ij,l,iq)
400c            print *,'q-->',ij,l,iq,q(ij,l,iq)
401             q(ij,l,iq)=zq(ij,l,iq)
402           ENDDO
403        ENDDO
404c$OMP END DO NOWAIT         
405
406c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
407        DO l=1,llm
408           DO ij=ijb,ije-iip1+1,iip1
409              q(ij+iim,l,iq)=q(ij,l,iq)
410           ENDDO
411        ENDDO
412c$OMP END DO NOWAIT 
413
414      ENDDO
415       
416     
417c$OMP BARRIER
418
419cc$OMP MASTER     
420c      call WaitSendRequest(MyRequest1)
421c      call WaitSendRequest(MyRequest2)
422cc$OMP END MASTER
423cc$OMP BARRIER
424
425      RETURN
426      END
Note: See TracBrowser for help on using the repository browser.