source: LMDZ5/branches/testing/libf/dyn3dpar/vlspltgen_p.F @ 4057

Last change on this file since 4057 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
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   --------------------------------------------------------------------
[1864]24      USE parallel_lmdz
[763]25      USE mod_hallo
26      USE Write_Field_p
27      USE VAMPIR
[1146]28      USE infotrac, ONLY : nqtot
[2641]29      USE comconst_mod, ONLY: cpp
[763]30      IMPLICIT NONE
31
32c
33#include "dimensions.h"
34#include "paramet.h"
35
36c
37c   Arguments:
38c   ----------
[1146]39      INTEGER iadv(nqtot)
[763]40      REAL masse(ip1jmp1,llm),pente_max
41      REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
[1146]42      REAL q(ip1jmp1,llm,nqtot)
[763]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)
[1146]52      REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: zm
[763]53      REAL,SAVE :: mu(ip1jmp1,llm)
54      REAL,SAVE :: mv(ip1jm,llm)
55      REAL,SAVE :: mw(ip1jmp1,llm+1)
[1146]56      REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: zq
[763]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
[1146]68      LOGICAL, SAVE :: firstcall=.TRUE.
69!$OMP THREADPRIVATE(firstcall)
[763]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
[1146]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
[763]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.
[985]100
[763]101      call SetTag(MyRequest1,100)
102      call SetTag(MyRequest2,101)
[985]103
[763]104       
[2641]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       
[763]110c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[2641]111        DO l = 1, llm
[763]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
[1146]177      DO iq=1,nqtot
[763]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           
[1146]188      DO iq=1,nqtot
[763]189
190        if(iadv(iq) == 0) then
[2641]191       
192          cycle
193       
194        else if (iadv(iq)==10) then
[985]195
196#ifdef _ADV_HALO       
197          call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
[2641]198     &                     ij_begin,ij_begin+2*iip1-1)
199          call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
[985]200     &               ij_end-2*iip1+1,ij_end)
201#else
[2641]202          call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
203     &                     ij_begin,ij_end)
[985]204#endif
[763]205
206c$OMP MASTER
207          call VTb(VTHallo)
[985]208c$OMP END MASTER
[763]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)
[985]211
212c$OMP MASTER
[763]213          call VTe(VTHallo)
214c$OMP END MASTER
[2641]215        else if (iadv(iq)==14) then
[985]216
217#ifdef _ADV_HALO           
[763]218          call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat,
[985]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,
[854]225     &                 ij_begin,ij_end)
[985]226#endif
227
[763]228c$OMP MASTER
229          call VTb(VTHallo)
[985]230c$OMP END MASTER
231
[763]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)
[985]234
235c$OMP MASTER
[763]236          call VTe(VTHallo)
237c$OMP END MASTER
238        else
[2641]239       
240          stop 'vlspltgen_p : schema non parallelise'
[763]241     
242        endif
243     
244      enddo
245     
246     
247c$OMP BARRIER     
248c$OMP MASTER     
249      call VTb(VTHallo)
[985]250c$OMP END MASTER
251
[763]252      call SendRequest(MyRequest1)
[985]253
254c$OMP MASTER
[763]255      call VTe(VTHallo)
[985]256c$OMP END MASTER       
257c$OMP BARRIER
[1146]258      do iq=1,nqtot
[985]259
260        if(iadv(iq) == 0) then
[2641]261       
262          cycle
263       
264        else if (iadv(iq)==10) then
[985]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       
[2641]270        else if (iadv(iq)==14) then
[985]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
[2641]276       
277          stop 'vlspltgen_p : schema non parallelise'
[985]278     
279        endif
280     
281      enddo
282c$OMP BARRIER     
283c$OMP MASTER
284      call VTb(VTHallo)
[763]285c$OMP END MASTER
[985]286
287!      call WaitRecvRequest(MyRequest1)
288!      call WaitSendRequest(MyRequest1)
[763]289c$OMP BARRIER
[985]290       call WaitRequest(MyRequest1)
291
292
293c$OMP MASTER
294      call VTe(VTHallo)
295c$OMP END MASTER
296c$OMP BARRIER
[763]297 
[1146]298      do iq=1,nqtot
[763]299
300        if(iadv(iq) == 0) then
301       
[2641]302          cycle
303       
304        else if (iadv(iq)==10) then
305       
[763]306          call vly_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv)
307 
[2641]308        else if (iadv(iq)==14) then
[763]309     
310          call vlyqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv,qsat)
311 
312        else
[2641]313       
314          stop 'vlspltgen_p : schema non parallelise'
[763]315     
316        endif
317       
318       enddo
319
[985]320
[1146]321      do iq=1,nqtot
[763]322
323        if(iadv(iq) == 0) then
[2641]324         
325          cycle
326       
327        else if (iadv(iq)==10 .or. iadv(iq)==14 ) then
[763]328
329c$OMP BARRIER       
[985]330#ifdef _ADV_HALLO
[763]331          call vlz_p(zq(1,1,iq),pente_max,zm(1,1,iq),mw,
[985]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,
[854]337     &               ij_begin,ij_end)
[985]338#endif
[763]339c$OMP BARRIER
340
341c$OMP MASTER
342          call VTb(VTHallo)
[985]343c$OMP END MASTER
344
[763]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)
[985]347
348c$OMP MASTER
[763]349          call VTe(VTHallo)
[2641]350c$OMP END MASTER       
[763]351c$OMP BARRIER
352        else
[2641]353       
354          stop 'vlspltgen_p : schema non parallelise'
[763]355     
356        endif
357     
358      enddo
359c$OMP BARRIER     
[985]360
[763]361c$OMP MASTER       
362      call VTb(VTHallo)
[985]363c$OMP END MASTER
364
[763]365      call SendRequest(MyRequest2)
[985]366
367c$OMP MASTER
[763]368      call VTe(VTHallo)
[2641]369c$OMP END MASTER       
[985]370
371c$OMP BARRIER
[1146]372      do iq=1,nqtot
[985]373
374        if(iadv(iq) == 0) then
[2641]375         
376          cycle
377       
378        else if (iadv(iq)==10 .or. iadv(iq)==14 ) then
[985]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
[2641]388       
389          stop 'vlspltgen_p : schema non parallelise'
[985]390     
391        endif
392     
393      enddo
394
395c$OMP BARRIER
396c$OMP MASTER
397      call VTb(VTHallo)
[763]398c$OMP END MASTER
[985]399
400!      call WaitRecvRequest(MyRequest2)
401!      call WaitSendRequest(MyRequest2)
[763]402c$OMP BARRIER
[985]403       CALL WaitRequest(MyRequest2)
[763]404
[985]405c$OMP MASTER
406      call VTe(VTHallo)
407c$OMP END MASTER
408c$OMP BARRIER
[763]409
[985]410
[1146]411      do iq=1,nqtot
[763]412
413        if(iadv(iq) == 0) then
414       
[2641]415          cycle
416       
417        else if (iadv(iq)==10) then
418       
[763]419          call vly_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv)
420 
[2641]421        else if (iadv(iq)==14) then
[763]422     
423          call vlyqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv,qsat)
424 
425        else
[2641]426       
427          stop 'vlspltgen_p : schema non parallelise'
[763]428     
429        endif
430       
431       enddo
432
[1146]433      do iq=1,nqtot
[763]434
435        if(iadv(iq) == 0) then
[2641]436         
437          cycle
[763]438       
[2641]439        else if (iadv(iq)==10) then
440       
[763]441          call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
442     &               ij_begin,ij_end)
443 
[2641]444        else if (iadv(iq)==14) then
[763]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
[2641]450       
[763]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
[985]462
[1146]463      DO iq=1,nqtot
[763]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)
[2641]469c             print *,'q-->',ij,l,iq,q(ij,l,iq)
470             q(ij,l,iq)=zq(ij,l,iq)
[763]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.