source: LMDZ5/branches/testing/libf/dyn3dmem/vlspltgen_loc.F @ 1790

Last change on this file since 1790 was 1669, checked in by Laurent Fairhead, 12 years ago

Version testing basée sur la r1668

http://lmdz.lmd.jussieu.fr/utilisateurs/distribution-du-modele/versions-intermediaires

Testing release based on r1668
File size: 11.5 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
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) :: 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 Allocate variables depending on dynamic variable nqtot
87
88         IF (firstcall) THEN
89            firstcall=.FALSE.
90         END IF
91c-- Calcul de Qsat en chaque point
92c-- approximation: au milieu des couches play(l)=(p(l)+p(l+1))/2
93c   pour eviter une exponentielle.
94
95      call SetTag(MyRequest1,100)
96      call SetTag(MyRequest2,101)
97
98       
99        ijb=ij_begin-iip1
100        ije=ij_end+iip1
101        if (pole_nord) ijb=ij_begin
102        if (pole_sud) ije=ij_end
103       
104c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
105        DO l = 1, llm
106         DO ij = ijb, ije
107          tempe(ij) = teta(ij,l) * pk(ij,l) /cpp
108         ENDDO
109         DO ij = ijb, ije
110          zdelta = MAX( 0., SIGN(1., rtt - tempe(ij)) )
111          play   = 0.5*(p(ij,l)+p(ij,l+1))
112          qsat(ij,l) = MIN(0.5, r2es* FOEEW(tempe(ij),zdelta) / play )
113          qsat(ij,l) = qsat(ij,l) / ( 1. - retv * qsat(ij,l) )
114         ENDDO
115        ENDDO
116c$OMP END DO NOWAIT
117c      PRINT*,'Debut vlsplt version debug sans vlyqs'
118
119        zzpbar = 0.5 * pdt
120        zzw    = pdt
121
122      ijb=ij_begin
123      ije=ij_end
124      if (pole_nord) ijb=ijb+iip1
125      if (pole_sud)  ije=ije-iip1
126
127c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
128      DO l=1,llm
129        DO ij = ijb,ije
130            mu(ij,l)=pbaru(ij,l) * zzpbar
131         ENDDO
132      ENDDO
133c$OMP END DO NOWAIT
134     
135      ijb=ij_begin-iip1
136      ije=ij_end
137      if (pole_nord) ijb=ij_begin
138      if (pole_sud)  ije=ij_end-iip1
139
140c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
141      DO l=1,llm
142         DO ij=ijb,ije
143            mv(ij,l)=pbarv(ij,l) * zzpbar
144         ENDDO
145      ENDDO
146c$OMP END DO NOWAIT
147
148      ijb=ij_begin
149      ije=ij_end
150
151c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
152      DO l=1,llm
153         DO ij=ijb,ije
154            mw(ij,l)=w(ij,l) * zzw
155         ENDDO
156      ENDDO
157c$OMP END DO NOWAIT
158
159c$OMP MASTER
160      DO ij=ijb,ije
161         mw(ij,llm+1)=0.
162      ENDDO
163c$OMP END MASTER
164
165c      CALL SCOPY(ijp1llm,q,1,zq,1)
166c      CALL SCOPY(ijp1llm,masse,1,zm,1)
167
168       ijb=ij_begin
169       ije=ij_end
170
171      DO iq=1,nqtot
172c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
173        DO l=1,llm
174          zq(ijb:ije,l,iq)=q(ijb:ije,l,iq)
175          zm(ijb:ije,l,iq)=masse(ijb:ije,l)
176        ENDDO
177c$OMP END DO NOWAIT
178      ENDDO
179
180#ifdef DEBUG_IO   
181       CALL WriteField_u('mu',mu)
182       CALL WriteField_v('mv',mv)
183       CALL WriteField_u('mw',mw)
184       CALL WriteField_u('qsat',qsat)
185#endif
186
187c$OMP BARRIER           
188      DO iq=1,nqtot
189
190#ifdef DEBUG_IO   
191       CALL WriteField_u('zq',zq(:,:,iq))
192       CALL WriteField_u('zm',zm(:,:,iq))
193#endif
194        if(iadv(iq) == 0) then
195       
196          cycle
197       
198        else if (iadv(iq)==10) then
199
200#ifdef _ADV_HALO       
201          call vlx_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
202     &               ij_begin,ij_begin+2*iip1-1)
203          call vlx_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
204     &               ij_end-2*iip1+1,ij_end)
205#else
206          call vlx_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
207     &               ij_begin,ij_end)
208#endif
209
210c$OMP MASTER
211          call VTb(VTHallo)
212c$OMP END MASTER
213          call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest1)
214          call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest1)
215
216c$OMP MASTER
217          call VTe(VTHallo)
218c$OMP END MASTER
219        else if (iadv(iq)==14) then
220
221#ifdef _ADV_HALO           
222          call vlxqs_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
223     &                   qsat,ij_begin,ij_begin+2*iip1-1)
224          call vlxqs_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
225     &                   qsat,ij_end-2*iip1+1,ij_end)
226#else
227
228          call vlxqs_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
229     &                   qsat,ij_begin,ij_end)
230#endif
231
232c$OMP MASTER
233          call VTb(VTHallo)
234c$OMP END MASTER
235
236          call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest1)
237          call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest1)
238
239c$OMP MASTER
240          call VTe(VTHallo)
241c$OMP END MASTER
242        else
243       
244          stop 'vlspltgen_p : schema non parallelise'
245     
246        endif
247     
248      enddo
249     
250     
251c$OMP BARRIER     
252c$OMP MASTER     
253      call VTb(VTHallo)
254c$OMP END MASTER
255
256      call SendRequest(MyRequest1)
257
258c$OMP MASTER
259      call VTe(VTHallo)
260c$OMP END MASTER       
261c$OMP BARRIER
262      do iq=1,nqtot
263
264        if(iadv(iq) == 0) then
265       
266          cycle
267       
268        else if (iadv(iq)==10) then
269
270#ifdef _ADV_HALLO
271          call vlx_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
272     &                 ij_begin+2*iip1,ij_end-2*iip1)
273#endif       
274        else if (iadv(iq)==14) then
275#ifdef _ADV_HALLO
276          call vlxqs_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
277     &                    qsat,ij_begin+2*iip1,ij_end-2*iip1)
278#endif   
279        else
280       
281          stop 'vlspltgen_p : schema non parallelise'
282     
283        endif
284     
285      enddo
286c$OMP BARRIER     
287c$OMP MASTER
288      call VTb(VTHallo)
289c$OMP END MASTER
290
291!      call WaitRecvRequest(MyRequest1)
292!      call WaitSendRequest(MyRequest1)
293c$OMP BARRIER
294       call WaitRequest(MyRequest1)
295
296
297c$OMP MASTER
298      call VTe(VTHallo)
299c$OMP END MASTER
300c$OMP BARRIER
301
302
303      do iq=1,nqtot
304#ifdef DEBUG_IO   
305       CALL WriteField_u('zq',zq(:,:,iq))
306       CALL WriteField_u('zm',zm(:,:,iq))
307#endif
308        if(iadv(iq) == 0) then
309       
310          cycle
311       
312        else if (iadv(iq)==10) then
313       
314          call vly_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mv)
315 
316        else if (iadv(iq)==14) then
317     
318          call vlyqs_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mv,
319     &                   qsat)
320 
321        else
322       
323          stop 'vlspltgen_p : schema non parallelise'
324     
325        endif
326       
327       enddo
328
329
330      do iq=1,nqtot
331#ifdef DEBUG_IO   
332       CALL WriteField_u('zq',zq(:,:,iq))
333       CALL WriteField_u('zm',zm(:,:,iq))
334#endif
335        if(iadv(iq) == 0) then
336         
337          cycle
338       
339        else if (iadv(iq)==10 .or. iadv(iq)==14 ) then
340
341c$OMP BARRIER       
342#ifdef _ADV_HALLO
343          call vlz_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mw,
344     &               ij_begin,ij_begin+2*iip1-1)
345          call vlz_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mw,
346     &               ij_end-2*iip1+1,ij_end)
347#else
348          call vlz_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mw,
349     &               ij_begin,ij_end)
350#endif
351c$OMP BARRIER
352
353c$OMP MASTER
354          call VTb(VTHallo)
355c$OMP END MASTER
356
357          call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest2)
358          call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest2)
359
360c$OMP MASTER
361          call VTe(VTHallo)
362c$OMP END MASTER       
363c$OMP BARRIER
364        else
365       
366          stop 'vlspltgen_p : schema non parallelise'
367     
368        endif
369     
370      enddo
371c$OMP BARRIER     
372
373c$OMP MASTER       
374      call VTb(VTHallo)
375c$OMP END MASTER
376
377      call SendRequest(MyRequest2)
378
379c$OMP MASTER
380      call VTe(VTHallo)
381c$OMP END MASTER       
382
383c$OMP BARRIER
384      do iq=1,nqtot
385
386        if(iadv(iq) == 0) then
387         
388          cycle
389       
390        else if (iadv(iq)==10 .or. iadv(iq)==14 ) then
391c$OMP BARRIER       
392
393#ifdef _ADV_HALLO
394          call vlz_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mw,
395     &               ij_begin+2*iip1,ij_end-2*iip1)
396#endif
397
398c$OMP BARRIER       
399        else
400       
401          stop 'vlspltgen_p : schema non parallelise'
402     
403        endif
404     
405      enddo
406
407c$OMP BARRIER
408c$OMP MASTER
409      call VTb(VTHallo)
410c$OMP END MASTER
411
412!      call WaitRecvRequest(MyRequest2)
413!      call WaitSendRequest(MyRequest2)
414c$OMP BARRIER
415       CALL WaitRequest(MyRequest2)
416
417c$OMP MASTER
418      call VTe(VTHallo)
419c$OMP END MASTER
420c$OMP BARRIER
421
422
423      do iq=1,nqtot
424#ifdef DEBUG_IO   
425       CALL WriteField_u('zq',zq(:,:,iq))
426       CALL WriteField_u('zm',zm(:,:,iq))
427#endif
428        if(iadv(iq) == 0) then
429       
430          cycle
431       
432        else if (iadv(iq)==10) then
433       
434          call vly_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mv)
435 
436        else if (iadv(iq)==14) then
437     
438          call vlyqs_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mv,
439     &                   qsat)
440 
441        else
442       
443          stop 'vlspltgen_p : schema non parallelise'
444     
445        endif
446       
447       enddo
448
449
450      do iq=1,nqtot
451#ifdef DEBUG_IO   
452       CALL WriteField_u('zq',zq(:,:,iq))
453       CALL WriteField_u('zm',zm(:,:,iq))
454#endif
455        if(iadv(iq) == 0) then
456         
457          cycle
458       
459        else if (iadv(iq)==10) then
460       
461          call vlx_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
462     &               ij_begin,ij_end)
463 
464        else if (iadv(iq)==14) then
465     
466          call vlxqs_loc(zq(ijb_u,1,iq),pente_max,zm(ijb_u,1,iq),mu,
467     &                 qsat, ij_begin,ij_end)
468 
469        else
470       
471          stop 'vlspltgen_p : schema non parallelise'
472     
473        endif
474       
475       enddo
476
477     
478      ijb=ij_begin
479      ije=ij_end
480c$OMP BARRIER     
481
482
483      DO iq=1,nqtot
484#ifdef DEBUG_IO   
485       CALL WriteField_u('zq',zq(:,:,iq))
486       CALL WriteField_u('zm',zm(:,:,iq))
487#endif
488c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
489        DO l=1,llm
490           DO ij=ijb,ije
491c             print *,'zq-->',ij,l,iq,zq(ij,l,iq)
492c            print *,'q-->',ij,l,iq,q(ij,l,iq)
493             q(ij,l,iq)=zq(ij,l,iq)
494           ENDDO
495        ENDDO
496c$OMP END DO NOWAIT         
497
498c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
499        DO l=1,llm
500           DO ij=ijb,ije-iip1+1,iip1
501              q(ij+iim,l,iq)=q(ij,l,iq)
502           ENDDO
503        ENDDO
504c$OMP END DO NOWAIT 
505
506      ENDDO
507       
508     
509c$OMP BARRIER
510
511cc$OMP MASTER     
512c      call WaitSendRequest(MyRequest1)
513c      call WaitSendRequest(MyRequest2)
514cc$OMP END MASTER
515cc$OMP BARRIER
516
517      RETURN
518      END
Note: See TracBrowser for help on using the repository browser.