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

Last change on this file since 923 was 854, checked in by lsce, 17 years ago

AC + JG + YM + ACo : - Corrected bugs found on 16 and 24 CPUs on ccrt scalar machine platine

  • Added use of lmdz communicator in inca model
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 8.2 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,
[854]186     &               ij_begin,ij_end)
[763]187
188c$OMP MASTER
189          call VTb(VTHallo)
190          call Register_Hallo(zq(1,1,iq),ip1jmp1,llm,2,2,2,2,MyRequest1)
191          call Register_Hallo(zm(1,1,iq),ip1jmp1,llm,1,1,1,1,MyRequest1)
192          call VTe(VTHallo)
193c$OMP END MASTER
194        else if (iadv(iq)==14) then
195     
196          call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat,
[854]197     &                 ij_begin,ij_end)
[763]198c$OMP MASTER
199          call VTb(VTHallo)
200          call Register_Hallo(zq(1,1,iq),ip1jmp1,llm,2,2,2,2,MyRequest1)
201          call Register_Hallo(zm(1,1,iq),ip1jmp1,llm,1,1,1,1,MyRequest1)
202          call VTe(VTHallo)
203c$OMP END MASTER
204        else
205       
206          stop 'vlspltgen_p : schema non parallelise'
207     
208        endif
209     
210      enddo
211     
212     
213c$OMP BARRIER     
214c$OMP MASTER     
215      call VTb(VTHallo)
216      call SendRequest(MyRequest1)
217      call WaitRecvRequest(MyRequest1)
218      call WaitSendRequest(MyRequest1)
219      call VTe(VTHallo)
220c$OMP END MASTER
221c$OMP BARRIER
222 
223      do iq=1,nqmx
224
225        if(iadv(iq) == 0) then
226       
227          cycle
228       
229        else if (iadv(iq)==10) then
230       
231          call vly_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv)
232 
233        else if (iadv(iq)==14) then
234     
235          call vlyqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv,qsat)
236 
237        else
238       
239          stop 'vlspltgen_p : schema non parallelise'
240     
241        endif
242       
243       enddo
244
245      do iq=1,nqmx
246
247        if(iadv(iq) == 0) then
248         
249          cycle
250       
251        else if (iadv(iq)==10 .or. iadv(iq)==14 ) then
252
253c$OMP BARRIER       
254          call vlz_p(zq(1,1,iq),pente_max,zm(1,1,iq),mw,
[854]255     &               ij_begin,ij_end)
[763]256c$OMP BARRIER
257
258c$OMP MASTER
259          call VTb(VTHallo)
260          call Register_Hallo(zq(1,1,iq),ip1jmp1,llm,2,2,2,2,MyRequest2)
261          call Register_Hallo(zm(1,1,iq),ip1jmp1,llm,1,1,1,1,MyRequest2)
262          call VTe(VTHallo)
263c$OMP END MASTER       
264c$OMP BARRIER
265        else
266       
267          stop 'vlspltgen_p : schema non parallelise'
268     
269        endif
270     
271      enddo
272c$OMP BARRIER     
273c$OMP MASTER       
274      call VTb(VTHallo)
275      call SendRequest(MyRequest2)
276      call WaitRecvRequest(MyRequest2)
277      call WaitSendRequest(MyRequest2)
278      call VTe(VTHallo)
279c$OMP END MASTER
280c$OMP BARRIER
281
282
283      do iq=1,nqmx
284
285        if(iadv(iq) == 0) then
286       
287          cycle
288       
289        else if (iadv(iq)==10) then
290       
291          call vly_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv)
292 
293        else if (iadv(iq)==14) then
294     
295          call vlyqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv,qsat)
296 
297        else
298       
299          stop 'vlspltgen_p : schema non parallelise'
300     
301        endif
302       
303       enddo
304
305      do iq=1,nqmx
306
307        if(iadv(iq) == 0) then
308         
309          cycle
310       
311        else if (iadv(iq)==10) then
312       
313          call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
314     &               ij_begin,ij_end)
315 
316        else if (iadv(iq)==14) then
317     
318          call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat,
319     &                 ij_begin,ij_end)
320 
321        else
322       
323          stop 'vlspltgen_p : schema non parallelise'
324     
325        endif
326       
327       enddo
328
329     
330      ijb=ij_begin
331      ije=ij_end
332c$OMP BARRIER     
333
334      DO iq=1,nqmx
335
336c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
337        DO l=1,llm
338           DO ij=ijb,ije
339c             print *,'zq-->',ij,l,iq,zq(ij,l,iq)
340c            print *,'q-->',ij,l,iq,q(ij,l,iq)
341             q(ij,l,iq)=zq(ij,l,iq)
342           ENDDO
343        ENDDO
344c$OMP END DO NOWAIT         
345
346c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
347        DO l=1,llm
348           DO ij=ijb,ije-iip1+1,iip1
349              q(ij+iim,l,iq)=q(ij,l,iq)
350           ENDDO
351        ENDDO
352c$OMP END DO NOWAIT 
353
354      ENDDO
355       
356     
357c$OMP BARRIER
358
359cc$OMP MASTER     
360c      call WaitSendRequest(MyRequest1)
361c      call WaitSendRequest(MyRequest2)
362cc$OMP END MASTER
363cc$OMP BARRIER
364
365      RETURN
366      END
Note: See TracBrowser for help on using the repository browser.