source: LMDZ5/branches/IPSLCM6.0.10/libf/dyn3dpar/vlspltgen_p.F @ 4778

Last change on this file since 4778 was 2641, checked in by Laurent Fairhead, 8 years ago

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