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

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

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