source: LMDZ5/branches/LMDZ6_rc0/libf/dyn3dmem/vlspltgen_loc.F @ 5394

Last change on this file since 5394 was 1910, checked in by Laurent Fairhead, 11 years ago

Merged trunk changes r1860:1909 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
File size: 11.6 KB
Line 
1!
2! $Header$
3!
4       SUBROUTINE vlspltgen_loc( q,iadv,pente_max,masse,w,pbaru,pbarv,
5     &                           pdt, p,pk,teta                 )
6     
7c
8c     Auteurs:   P.Le Van, F.Hourdin, F.Forget, F.Codron
9c
10c    ********************************************************************
11c          Shema  d'advection " pseudo amont " .
12c      + test sur humidite specifique: Q advecte< Qsat aval
13c                   (F. Codron, 10/99)
14c    ********************************************************************
15c     q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
16c
17c     pente_max facteur de limitation des pentes: 2 en general
18c                                                0 pour un schema amont
19c     pbaru,pbarv,w flux de masse en u ,v ,w
20c     pdt pas de temps
21c
22c     teta temperature potentielle, p pression aux interfaces,
23c     pk exner au milieu des couches necessaire pour calculer Qsat
24c   --------------------------------------------------------------------
25      USE parallel_lmdz
26      USE mod_hallo
27      USE Write_Field_loc
28      USE VAMPIR
29      USE infotrac, ONLY : nqtot
30      USE vlspltgen_mod
31      IMPLICIT NONE
32
33c
34#include "dimensions.h"
35#include "paramet.h"
36#include "logic.h"
37#include "comvert.h"
38#include "comconst.h"
39
40c
41c   Arguments:
42c   ----------
43      INTEGER iadv(nqtot)
44      REAL masse(ijb_u:ije_u,llm),pente_max
45      REAL pbaru( ijb_u:ije_u,llm ),pbarv( ijb_v:ije_v,llm)
46      REAL q(ijb_u:ije_u,llm,nqtot)
47      REAL w(ijb_u:ije_u,llm),pdt
48      REAL p(ijb_u:ije_u,llmp1),teta(ijb_u:ije_u,llm)
49      REAL pk(ijb_u:ije_u,llm)
50c
51c      Local
52c   ---------
53c
54      INTEGER ij,l
55c
56      REAL zzpbar, zzw
57
58      REAL qmin,qmax
59      DATA qmin,qmax/0.,1.e33/
60
61c--pour rapport de melange saturant--
62
63      REAL rtt,retv,r2es,r3les,r3ies,r4les,r4ies,play
64      REAL ptarg,pdelarg,foeew,zdelta
65      REAL tempe(ijb_u:ije_u)
66      INTEGER ijb,ije,iq
67      LOGICAL, SAVE :: firstcall=.TRUE.
68!$OMP THREADPRIVATE(firstcall)
69      type(request),SAVE :: MyRequest1
70!$OMP THREADPRIVATE(MyRequest1)
71      type(request),SAVE :: MyRequest2
72!$OMP THREADPRIVATE(MyRequest2)
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         END IF
92c-- Calcul de Qsat en chaque point
93c-- approximation: au milieu des couches play(l)=(p(l)+p(l+1))/2
94c   pour eviter une exponentielle.
95
96      call SetTag(MyRequest1,100)
97      call SetTag(MyRequest2,101)
98
99       
100        ijb=ij_begin-iip1
101        ije=ij_end+iip1
102        if (pole_nord) ijb=ij_begin
103        if (pole_sud) ije=ij_end
104       
105c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
106        DO l = 1, llm
107         DO ij = ijb, ije
108          tempe(ij) = teta(ij,l) * pk(ij,l) /cpp
109         ENDDO
110         DO ij = ijb, ije
111          zdelta = MAX( 0., SIGN(1., rtt - tempe(ij)) )
112          play   = 0.5*(p(ij,l)+p(ij,l+1))
113          qsat(ij,l) = MIN(0.5, r2es* FOEEW(tempe(ij),zdelta) / play )
114          qsat(ij,l) = qsat(ij,l) / ( 1. - retv * qsat(ij,l) )
115         ENDDO
116        ENDDO
117c$OMP END DO NOWAIT
118c      PRINT*,'Debut vlsplt version debug sans vlyqs'
119
120        zzpbar = 0.5 * pdt
121        zzw    = pdt
122
123      ijb=ij_begin
124      ije=ij_end
125      if (pole_nord) ijb=ijb+iip1
126      if (pole_sud)  ije=ije-iip1
127
128c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
129      DO l=1,llm
130        DO ij = ijb,ije
131            mu(ij,l)=pbaru(ij,l) * zzpbar
132         ENDDO
133      ENDDO
134c$OMP END DO NOWAIT
135     
136      ijb=ij_begin-iip1
137      ije=ij_end
138      if (pole_nord) ijb=ij_begin
139      if (pole_sud)  ije=ij_end-iip1
140
141c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
142      DO l=1,llm
143         DO ij=ijb,ije
144            mv(ij,l)=pbarv(ij,l) * zzpbar
145         ENDDO
146      ENDDO
147c$OMP END DO NOWAIT
148
149      ijb=ij_begin
150      ije=ij_end
151
152c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
153      DO l=1,llm
154         DO ij=ijb,ije
155            mw(ij,l)=w(ij,l) * zzw
156         ENDDO
157      ENDDO
158c$OMP END DO NOWAIT
159
160c$OMP MASTER
161      DO ij=ijb,ije
162         mw(ij,llm+1)=0.
163      ENDDO
164c$OMP END MASTER
165
166c      CALL SCOPY(ijp1llm,q,1,zq,1)
167c      CALL SCOPY(ijp1llm,masse,1,zm,1)
168
169       ijb=ij_begin
170       ije=ij_end
171
172      DO iq=1,nqtot
173c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
174        DO l=1,llm
175          zq(ijb:ije,l,iq)=q(ijb:ije,l,iq)
176          zm(ijb:ije,l,iq)=masse(ijb:ije,l)
177        ENDDO
178c$OMP END DO NOWAIT
179      ENDDO
180
181#ifdef DEBUG_IO   
182       CALL WriteField_u('mu',mu)
183       CALL WriteField_v('mv',mv)
184       CALL WriteField_u('mw',mw)
185       CALL WriteField_u('qsat',qsat)
186#endif
187
188c$OMP BARRIER           
189      DO iq=1,nqtot
190
191#ifdef DEBUG_IO   
192       CALL WriteField_u('zq',zq(:,:,iq))
193       CALL WriteField_u('zm',zm(:,:,iq))
194#endif
195        if(iadv(iq) == 0) then
196       
197          cycle
198       
199        else if (iadv(iq)==10) then
200
201#ifdef _ADV_HALO       
202          call vlx_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
203     &               ij_begin,ij_begin+2*iip1-1)
204          call vlx_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
205     &               ij_end-2*iip1+1,ij_end)
206#else
207          call vlx_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
208     &               ij_begin,ij_end)
209#endif
210
211c$OMP MASTER
212          call VTb(VTHallo)
213c$OMP END MASTER
214          call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest1)
215          call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest1)
216
217c$OMP MASTER
218          call VTe(VTHallo)
219c$OMP END MASTER
220        else if (iadv(iq)==14) then
221
222#ifdef _ADV_HALO           
223          call vlxqs_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
224     &                   qsat,ij_begin,ij_begin+2*iip1-1)
225          call vlxqs_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
226     &                   qsat,ij_end-2*iip1+1,ij_end)
227#else
228
229          call vlxqs_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
230     &                   qsat,ij_begin,ij_end)
231#endif
232
233c$OMP MASTER
234          call VTb(VTHallo)
235c$OMP END MASTER
236
237          call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest1)
238          call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest1)
239
240c$OMP MASTER
241          call VTe(VTHallo)
242c$OMP END MASTER
243        else
244       
245          stop 'vlspltgen_p : schema non parallelise'
246     
247        endif
248     
249      enddo
250     
251     
252c$OMP BARRIER     
253c$OMP MASTER     
254      call VTb(VTHallo)
255c$OMP END MASTER
256
257      call SendRequest(MyRequest1)
258
259c$OMP MASTER
260      call VTe(VTHallo)
261c$OMP END MASTER       
262c$OMP BARRIER
263      do iq=1,nqtot
264
265        if(iadv(iq) == 0) then
266       
267          cycle
268       
269        else if (iadv(iq)==10) then
270
271#ifdef _ADV_HALLO
272          call vlx_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
273     &                 ij_begin+2*iip1,ij_end-2*iip1)
274#endif       
275        else if (iadv(iq)==14) then
276#ifdef _ADV_HALLO
277          call vlxqs_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
278     &                    qsat,ij_begin+2*iip1,ij_end-2*iip1)
279#endif   
280        else
281       
282          stop 'vlspltgen_p : schema non parallelise'
283     
284        endif
285     
286      enddo
287c$OMP BARRIER     
288c$OMP MASTER
289      call VTb(VTHallo)
290c$OMP END MASTER
291
292!      call WaitRecvRequest(MyRequest1)
293!      call WaitSendRequest(MyRequest1)
294c$OMP BARRIER
295       call WaitRequest(MyRequest1)
296
297
298c$OMP MASTER
299      call VTe(VTHallo)
300c$OMP END MASTER
301c$OMP BARRIER
302
303
304      do iq=1,nqtot
305#ifdef DEBUG_IO   
306       CALL WriteField_u('zq',zq(:,:,iq))
307       CALL WriteField_u('zm',zm(:,:,iq))
308#endif
309        if(iadv(iq) == 0) then
310       
311          cycle
312       
313        else if (iadv(iq)==10) then
314       
315          call vly_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mv)
316 
317        else if (iadv(iq)==14) then
318     
319          call vlyqs_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mv,
320     &                   qsat)
321 
322        else
323       
324          stop 'vlspltgen_p : schema non parallelise'
325     
326        endif
327       
328       enddo
329
330
331      do iq=1,nqtot
332#ifdef DEBUG_IO   
333       CALL WriteField_u('zq',zq(:,:,iq))
334       CALL WriteField_u('zm',zm(:,:,iq))
335#endif
336        if(iadv(iq) == 0) then
337         
338          cycle
339       
340        else if (iadv(iq)==10 .or. iadv(iq)==14 ) then
341
342c$OMP BARRIER       
343#ifdef _ADV_HALLO
344          call vlz_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mw,
345     &               ij_begin,ij_begin+2*iip1-1)
346          call vlz_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mw,
347     &               ij_end-2*iip1+1,ij_end)
348#else
349          call vlz_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mw,
350     &               ij_begin,ij_end)
351#endif
352c$OMP BARRIER
353
354c$OMP MASTER
355          call VTb(VTHallo)
356c$OMP END MASTER
357
358          call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest2)
359          call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest2)
360
361c$OMP MASTER
362          call VTe(VTHallo)
363c$OMP END MASTER       
364c$OMP BARRIER
365        else
366       
367          stop 'vlspltgen_p : schema non parallelise'
368     
369        endif
370     
371      enddo
372c$OMP BARRIER     
373
374c$OMP MASTER       
375      call VTb(VTHallo)
376c$OMP END MASTER
377
378      call SendRequest(MyRequest2)
379
380c$OMP MASTER
381      call VTe(VTHallo)
382c$OMP END MASTER       
383
384c$OMP BARRIER
385      do iq=1,nqtot
386
387        if(iadv(iq) == 0) then
388         
389          cycle
390       
391        else if (iadv(iq)==10 .or. iadv(iq)==14 ) then
392c$OMP BARRIER       
393
394#ifdef _ADV_HALLO
395          call vlz_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mw,
396     &               ij_begin+2*iip1,ij_end-2*iip1)
397#endif
398
399c$OMP BARRIER       
400        else
401       
402          stop 'vlspltgen_p : schema non parallelise'
403     
404        endif
405     
406      enddo
407
408c$OMP BARRIER
409c$OMP MASTER
410      call VTb(VTHallo)
411c$OMP END MASTER
412
413!      call WaitRecvRequest(MyRequest2)
414!      call WaitSendRequest(MyRequest2)
415c$OMP BARRIER
416       CALL WaitRequest(MyRequest2)
417
418c$OMP MASTER
419      call VTe(VTHallo)
420c$OMP END MASTER
421c$OMP BARRIER
422
423
424      do iq=1,nqtot
425#ifdef DEBUG_IO   
426       CALL WriteField_u('zq',zq(:,:,iq))
427       CALL WriteField_u('zm',zm(:,:,iq))
428#endif
429        if(iadv(iq) == 0) then
430       
431          cycle
432       
433        else if (iadv(iq)==10) then
434       
435          call vly_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mv)
436 
437        else if (iadv(iq)==14) then
438     
439          call vlyqs_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mv,
440     &                   qsat)
441 
442        else
443       
444          stop 'vlspltgen_p : schema non parallelise'
445     
446        endif
447       
448       enddo
449
450
451      do iq=1,nqtot
452#ifdef DEBUG_IO   
453       CALL WriteField_u('zq',zq(:,:,iq))
454       CALL WriteField_u('zm',zm(:,:,iq))
455#endif
456        if(iadv(iq) == 0) then
457         
458          cycle
459       
460        else if (iadv(iq)==10) then
461       
462          call vlx_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
463     &               ij_begin,ij_end)
464 
465        else if (iadv(iq)==14) then
466     
467          call vlxqs_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
468     &                 qsat, ij_begin,ij_end)
469 
470        else
471       
472          stop 'vlspltgen_p : schema non parallelise'
473     
474        endif
475       
476       enddo
477
478     
479      ijb=ij_begin
480      ije=ij_end
481c$OMP BARRIER     
482
483
484      DO iq=1,nqtot
485#ifdef DEBUG_IO   
486       CALL WriteField_u('zq',zq(:,:,iq))
487       CALL WriteField_u('zm',zm(:,:,iq))
488#endif
489c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
490        DO l=1,llm
491           DO ij=ijb,ije
492c             print *,'zq-->',ij,l,iq,zq(ij,l,iq)
493c            print *,'q-->',ij,l,iq,q(ij,l,iq)
494             q(ij,l,iq)=zq(ij,l,iq)
495           ENDDO
496        ENDDO
497c$OMP END DO NOWAIT         
498
499c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
500        DO l=1,llm
501           DO ij=ijb,ije-iip1+1,iip1
502              q(ij+iim,l,iq)=q(ij,l,iq)
503           ENDDO
504        ENDDO
505c$OMP END DO NOWAIT 
506
507      ENDDO
508       
509     
510c$OMP BARRIER
511
512cc$OMP MASTER     
513c      call WaitSendRequest(MyRequest1)
514c      call WaitSendRequest(MyRequest2)
515cc$OMP END MASTER
516cc$OMP BARRIER
517
518      RETURN
519      END
Note: See TracBrowser for help on using the repository browser.