source: trunk/LMDZ.COMMON/libf/dyn3dpar/vlspltgen_p.F @ 1685

Last change on this file since 1685 was 1422, checked in by milmd, 10 years ago

In GENERIC, MARS and COMMON models replace some include files by modules (usefull for decoupling physics with dynamics).
File size: 10.7 KB
RevLine 
[1]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   --------------------------------------------------------------------
[1019]24      USE parallel_lmdz
[1]25      USE mod_hallo
26      USE Write_Field_p
27      USE VAMPIR
28      USE infotrac, ONLY : nqtot
[1422]29      USE comconst_mod, ONLY: cpp
[1]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.