source: LMDZ5/trunk/libf/dyn3dmem/vlspltgen_loc.F @ 2278

Last change on this file since 2278 was 2270, checked in by crisi, 9 years ago

Adding isotopes in the dynamics and more generally tracers of tracers.
CRisi
  • 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: 14.9 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      ! CRisi: on rajoute variables utiles d'infotrac 
30      USE infotrac, ONLY : nqtot,nqperes,nqdesc,nqfils,iqfils,
31     &    ok_iso_verif
32      USE vlspltgen_mod
33      IMPLICIT NONE
34
35c
36#include "dimensions.h"
37#include "paramet.h"
38#include "logic.h"
39#include "comvert.h"
40#include "comconst.h"
41
42c
43c   Arguments:
44c   ----------
45      INTEGER iadv(nqtot)
46      REAL masse(ijb_u:ije_u,llm),pente_max
47      REAL pbaru( ijb_u:ije_u,llm ),pbarv( ijb_v:ije_v,llm)
48      REAL q(ijb_u:ije_u,llm,nqtot)
49      REAL w(ijb_u:ije_u,llm),pdt
50      REAL p(ijb_u:ije_u,llmp1),teta(ijb_u:ije_u,llm)
51      REAL pk(ijb_u:ije_u,llm)
52c
53c      Local
54c   ---------
55c
56      INTEGER ij,l
57c
58      REAL zzpbar, zzw
59
60      REAL qmin,qmax
61      DATA qmin,qmax/0.,1.e33/
62
63c--pour rapport de melange saturant--
64
65      REAL rtt,retv,r2es,r3les,r3ies,r4les,r4ies,play
66      REAL ptarg,pdelarg,foeew,zdelta
67      REAL tempe(ijb_u:ije_u)
68      INTEGER ijb,ije,iq,iq2,ifils
69      LOGICAL, SAVE :: firstcall=.TRUE.
70!$OMP THREADPRIVATE(firstcall)
71      type(request),SAVE :: MyRequest1
72!$OMP THREADPRIVATE(MyRequest1)
73      type(request),SAVE :: MyRequest2
74!$OMP THREADPRIVATE(MyRequest2)
75c    fonction psat(T)
76
77       FOEEW ( PTARG,PDELARG ) = EXP (
78     *          (R3LES*(1.-PDELARG)+R3IES*PDELARG) * (PTARG-RTT)
79     * / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG)) )
80
81        r2es  = 380.11733
82        r3les = 17.269
83        r3ies = 21.875
84        r4les = 35.86
85        r4ies = 7.66
86        retv = 0.6077667
87        rtt  = 273.16
88
89c Allocate variables depending on dynamic variable nqtot
90
91         IF (firstcall) THEN
92            firstcall=.FALSE.
93         END IF
94c-- Calcul de Qsat en chaque point
95c-- approximation: au milieu des couches play(l)=(p(l)+p(l+1))/2
96c   pour eviter une exponentielle.
97
98      call SetTag(MyRequest1,100)
99      call SetTag(MyRequest2,101)
100
101       
102        ijb=ij_begin-iip1
103        ije=ij_end+iip1
104        if (pole_nord) ijb=ij_begin
105        if (pole_sud) ije=ij_end
106       
107c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
108        DO l = 1, llm
109         DO ij = ijb, ije
110          tempe(ij) = teta(ij,l) * pk(ij,l) /cpp
111         ENDDO
112         DO ij = ijb, ije
113          zdelta = MAX( 0., SIGN(1., rtt - tempe(ij)) )
114          play   = 0.5*(p(ij,l)+p(ij,l+1))
115          qsat(ij,l) = MIN(0.5, r2es* FOEEW(tempe(ij),zdelta) / play )
116          qsat(ij,l) = qsat(ij,l) / ( 1. - retv * qsat(ij,l) )
117         ENDDO
118        ENDDO
119c$OMP END DO NOWAIT
120c      PRINT*,'Debut vlsplt version debug sans vlyqs'
121
122        zzpbar = 0.5 * pdt
123        zzw    = pdt
124
125      ijb=ij_begin
126      ije=ij_end
127      if (pole_nord) ijb=ijb+iip1
128      if (pole_sud)  ije=ije-iip1
129
130c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
131      DO l=1,llm
132        DO ij = ijb,ije
133            mu(ij,l)=pbaru(ij,l) * zzpbar
134         ENDDO
135      ENDDO
136c$OMP END DO NOWAIT
137     
138      ijb=ij_begin-iip1
139      ije=ij_end
140      if (pole_nord) ijb=ij_begin
141      if (pole_sud)  ije=ij_end-iip1
142
143c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
144      DO l=1,llm
145         DO ij=ijb,ije
146            mv(ij,l)=pbarv(ij,l) * zzpbar
147         ENDDO
148      ENDDO
149c$OMP END DO NOWAIT
150
151      ijb=ij_begin
152      ije=ij_end
153
154      DO iq=1,nqtot
155c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
156      DO l=1,llm
157         DO ij=ijb,ije
158            mw(ij,l,iq)=w(ij,l) * zzw
159         ENDDO
160      ENDDO
161c$OMP END DO NOWAIT
162      ENDDO
163
164      DO iq=1,nqtot 
165c$OMP MASTER
166      DO ij=ijb,ije
167         mw(ij,llm+1,iq)=0.
168      ENDDO
169c$OMP END MASTER
170      ENDDO
171
172c      CALL SCOPY(ijp1llm,q,1,zq,1)
173c      CALL SCOPY(ijp1llm,masse,1,zm,1)
174
175       ijb=ij_begin
176       ije=ij_end
177
178      DO iq=1,nqtot       
179c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
180        DO l=1,llm
181          zq(ijb:ije,l,iq)=q(ijb:ije,l,iq)
182          zm(ijb:ije,l,iq)=masse(ijb:ije,l)
183        ENDDO
184c$OMP END DO NOWAIT
185      ENDDO
186
187#ifdef DEBUG_IO     
188       CALL WriteField_u('mu',mu)
189       CALL WriteField_v('mv',mv)
190       CALL WriteField_u('mw',mw)
191       CALL WriteField_u('qsat',qsat)
192#endif
193
194      ! verif temporaire
195      ijb=ij_begin-2*iip1
196      ije=ij_end+2*iip1 
197      if (pole_nord) ijb=ij_begin
198      if (pole_sud)  ije=ij_end 
199      if (ok_iso_verif) then
200           call check_isotopes(zq,ijb,ije,'vlspltgen_loc 191')
201      endif !if (ok_iso_verif) then   
202
203c$OMP BARRIER           
204!      DO iq=1,nqtot
205      DO iq=1,nqperes ! CRisi: on ne boucle que sur les pères= ceux qui sont transportés directement par l'air
206       write(*,*) 'vlspltgen 192: iq,iadv=',iq,iadv(iq)
207#ifdef DEBUG_IO   
208       CALL WriteField_u('zq',zq(:,:,iq))
209       CALL WriteField_u('zm',zm(:,:,iq))
210#endif
211        if(iadv(iq) == 0) then
212       
213          cycle
214       
215        else if (iadv(iq)==10) then
216
217#ifdef _ADV_HALO       
218! CRisi: on ajoute les nombres de fils et tableaux des fils
219! On suppose qu'on ne peut advecter les fils que par le schéma 10. 
220          call vlx_loc(zq,pente_max,zm,mu,
221     &               ij_begin,ij_begin+2*iip1-1,iq)
222          call vlx_loc(zq,pente_max,zm,mu,
223     &               ij_end-2*iip1+1,ij_end,iq)
224#else
225          call vlx_loc(zq,pente_max,zm,mu,
226     &               ij_begin,ij_end,iq)
227#endif
228
229c$OMP MASTER
230          call VTb(VTHallo)
231c$OMP END MASTER
232          call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest1)
233          call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest1)
234! CRisi
235          do ifils=1,nqdesc(iq)
236            iq2=iqfils(ifils,iq)
237            call Register_Hallo_u(zq(:,:,iq2),llm,2,2,2,2,MyRequest1)
238            call Register_Hallo_u(zm(:,:,iq2),llm,1,1,1,1,MyRequest1)
239          enddo
240
241c$OMP MASTER
242          call VTe(VTHallo)
243c$OMP END MASTER
244        else if (iadv(iq)==14) then
245
246#ifdef _ADV_HALO           
247          call vlxqs_loc(zq,pente_max,zm,mu,
248     &                   qsat,ij_begin,ij_begin+2*iip1-1,iq)
249          call vlxqs_loc(zq,pente_max,zm,mu,
250     &                   qsat,ij_end-2*iip1+1,ij_end,iq)
251#else
252          call vlxqs_loc(zq,pente_max,zm,mu,
253     &                   qsat,ij_begin,ij_end,iq)
254#endif
255
256c$OMP MASTER
257          call VTb(VTHallo)
258c$OMP END MASTER
259
260          call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest1)
261          call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest1)
262          do ifils=1,nqdesc(iq)
263            iq2=iqfils(ifils,iq)
264            call Register_Hallo_u(zq(:,:,iq2),llm,2,2,2,2,MyRequest1)
265            call Register_Hallo_u(zm(:,:,iq2),llm,1,1,1,1,MyRequest1)
266          enddo
267
268c$OMP MASTER
269          call VTe(VTHallo)
270c$OMP END MASTER
271        else
272       
273          stop 'vlspltgen_p : schema non parallelise'
274     
275        endif
276     
277      enddo !DO iq=1,nqperes
278     
279     
280c$OMP BARRIER     
281c$OMP MASTER     
282      call VTb(VTHallo)
283c$OMP END MASTER
284
285      call SendRequest(MyRequest1)
286
287c$OMP MASTER
288      call VTe(VTHallo)
289c$OMP END MASTER       
290c$OMP BARRIER
291
292      ! verif temporaire
293      ijb=ij_begin-2*iip1
294      ije=ij_end+2*iip1 
295      if (pole_nord) ijb=ij_begin
296      if (pole_sud)  ije=ij_end 
297      if (ok_iso_verif) then
298           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 280')
299      endif !if (ok_iso_verif) then
300
301      do iq=1,nqperes
302        write(*,*) 'vlspltgen 279: iq=',iq
303
304        if(iadv(iq) == 0) then
305       
306          cycle
307       
308        else if (iadv(iq)==10) then
309
310#ifdef _ADV_HALLO
311          call vlx_loc(zq,pente_max,zm,mu,
312     &                 ij_begin+2*iip1,ij_end-2*iip1,iq)
313#endif       
314        else if (iadv(iq)==14) then
315#ifdef _ADV_HALLO
316          call vlxqs_loc(zq,pente_max,zm,mu,
317     &                    qsat,ij_begin+2*iip1,ij_end-2*iip1,iq)
318#endif   
319        else
320       
321          stop 'vlspltgen_p : schema non parallelise'
322     
323        endif
324     
325      enddo
326c$OMP BARRIER     
327c$OMP MASTER
328      call VTb(VTHallo)
329c$OMP END MASTER
330
331!      call WaitRecvRequest(MyRequest1)
332!      call WaitSendRequest(MyRequest1)
333c$OMP BARRIER
334       call WaitRequest(MyRequest1)
335
336
337c$OMP MASTER
338      call VTe(VTHallo)
339c$OMP END MASTER
340c$OMP BARRIER
341
342     
343      if (ok_iso_verif) then
344           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 326')
345      endif !if (ok_iso_verif) then       
346      if (ok_iso_verif) then
347           ijb=ij_begin-2*iip1
348           ije=ij_end+2*iip1
349           if (pole_nord) ijb=ij_begin
350           if (pole_sud)  ije=ij_end
351           call check_isotopes(zq,ijb,ije,'vlspltgen_loc 336')
352      endif !if (ok_iso_verif) then 
353
354      do iq=1,nqperes
355       write(*,*) 'vlspltgen 321: iq=',iq
356#ifdef DEBUG_IO   
357       CALL WriteField_u('zq',zq(:,:,iq))
358       CALL WriteField_u('zm',zm(:,:,iq))
359#endif
360
361        if(iadv(iq) == 0) then
362       
363          cycle
364       
365        else if (iadv(iq)==10) then
366       
367          call vly_loc(zq,pente_max,zm,mv,iq)
368 
369        else if (iadv(iq)==14) then
370     
371          call vlyqs_loc(zq,pente_max,zm,mv,
372     &                   qsat,iq)
373 
374        else
375       
376          stop 'vlspltgen_p : schema non parallelise'
377     
378        endif
379       
380       enddo
381
382      if (ok_iso_verif) then
383           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 357')
384      endif !if (ok_iso_verif) then
385
386      do iq=1,nqperes
387      write(*,*) 'vlspltgen 349: iq=',iq
388#ifdef DEBUG_IO   
389       CALL WriteField_u('zq',zq(:,:,iq))
390       CALL WriteField_u('zm',zm(:,:,iq))
391#endif
392        if(iadv(iq) == 0) then
393         
394          cycle
395       
396        else if (iadv(iq)==10 .or. iadv(iq)==14 ) then
397
398c$OMP BARRIER       
399#ifdef _ADV_HALLO
400          call vlz_loc(zq,pente_max,zm,mw,
401     &               ij_begin,ij_begin+2*iip1-1,iq)
402          call vlz_loc(zq,pente_max,zm,mw,
403     &               ij_end-2*iip1+1,ij_end,iq)
404#else
405          call vlz_loc(zq,pente_max,zm,mw,
406     &               ij_begin,ij_end,iq)
407#endif
408c$OMP BARRIER
409
410c$OMP MASTER
411          call VTb(VTHallo)
412c$OMP END MASTER
413
414          call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest2)
415          call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest2)
416          ! CRisi
417          do ifils=1,nqdesc(iq)
418            iq2=iqfils(ifils,iq)
419            call Register_Hallo_u(zq(:,:,iq2),llm,2,2,2,2,MyRequest2)
420            call Register_Hallo_u(zm(:,:,iq2),llm,1,1,1,1,MyRequest2)
421          enddo     
422c$OMP MASTER
423          call VTe(VTHallo)
424c$OMP END MASTER       
425c$OMP BARRIER
426        else
427       
428          stop 'vlspltgen_p : schema non parallelise'
429     
430        endif
431     
432      enddo
433c$OMP BARRIER     
434
435c$OMP MASTER       
436      call VTb(VTHallo)
437c$OMP END MASTER
438
439      call SendRequest(MyRequest2)
440
441c$OMP MASTER
442      call VTe(VTHallo)
443c$OMP END MASTER       
444
445
446      if (ok_iso_verif) then
447           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 429')
448      endif !if (ok_iso_verif) then
449
450c$OMP BARRIER
451      do iq=1,nqperes
452      write(*,*) 'vlspltgen 409: iq=',iq
453
454        if(iadv(iq) == 0) then
455         
456          cycle
457       
458        else if (iadv(iq)==10 .or. iadv(iq)==14 ) then
459c$OMP BARRIER       
460
461          write(*,*) 'vlspltgen_loc 461'
462#ifdef _ADV_HALLO
463          write(*,*) 'vlspltgen_loc 462'
464          call vlz_loc(zq,pente_max,zm,mw,
465     &               ij_begin+2*iip1,ij_end-2*iip1,iq)
466#endif
467
468c$OMP BARRIER       
469        else
470       
471          stop 'vlspltgen_p : schema non parallelise'
472     
473        endif
474     
475      enddo
476      write(*,*) 'vlspltgen_loc 476'
477
478c$OMP BARRIER
479c$OMP MASTER
480      call VTb(VTHallo)
481c$OMP END MASTER
482
483!      call WaitRecvRequest(MyRequest2)
484!      call WaitSendRequest(MyRequest2)
485c$OMP BARRIER
486       CALL WaitRequest(MyRequest2)
487
488c$OMP MASTER
489      call VTe(VTHallo)
490c$OMP END MASTER
491c$OMP BARRIER
492
493
494      write(*,*) 'vlspltgen_loc 494'
495      if (ok_iso_verif) then
496           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 461')
497      endif !if (ok_iso_verif) then
498
499      do iq=1,nqperes
500      write(*,*) 'vlspltgen 449: iq=',iq
501#ifdef DEBUG_IO   
502       CALL WriteField_u('zq',zq(:,:,iq))
503       CALL WriteField_u('zm',zm(:,:,iq))
504#endif
505        if(iadv(iq) == 0) then
506       
507          cycle
508       
509        else if (iadv(iq)==10) then
510       
511          call vly_loc(zq,pente_max,zm,mv,iq)
512 
513        else if (iadv(iq)==14) then
514     
515          call vlyqs_loc(zq,pente_max,zm,mv,
516     &                   qsat,iq)
517 
518        else
519       
520          stop 'vlspltgen_p : schema non parallelise'
521     
522        endif
523       
524       enddo !do iq=1,nqperes
525
526      if (ok_iso_verif) then
527           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 493')
528      endif !if (ok_iso_verif) then
529
530      do iq=1,nqperes
531      write(*,*) 'vlspltgen 477: iq=',iq
532#ifdef DEBUG_IO   
533       CALL WriteField_u('zq',zq(:,:,iq))
534       CALL WriteField_u('zm',zm(:,:,iq))
535#endif
536        if(iadv(iq) == 0) then
537         
538          cycle
539       
540        else if (iadv(iq)==10) then
541       
542          call vlx_loc(zq,pente_max,zm,mu,
543     &               ij_begin,ij_end,iq)
544 
545        else if (iadv(iq)==14) then
546     
547          call vlxqs_loc(zq,pente_max,zm,mu,
548     &                 qsat, ij_begin,ij_end,iq)
549 
550        else
551       
552          stop 'vlspltgen_p : schema non parallelise'
553     
554        endif
555       
556       enddo !do iq=1,nqperes
557
558      write(*,*) 'vlspltgen 550: apres derniere serie de call vlx'
559      if (ok_iso_verif) then
560           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 521')
561      endif !if (ok_iso_verif) then
562     
563      ijb=ij_begin
564      ije=ij_end
565      write(*,*) 'vlspltgen_loc 557'
566c$OMP BARRIER     
567
568      write(*,*) 'vlspltgen_loc 559' 
569      DO iq=1,nqtot
570       write(*,*) 'vlspltgen_loc 561, iq=',iq 
571#ifdef DEBUG_IO   
572       CALL WriteField_u('zq',zq(:,:,iq))
573       CALL WriteField_u('zm',zm(:,:,iq))
574#endif
575c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
576        DO l=1,llm
577           DO ij=ijb,ije
578c             print *,'zq-->',ij,l,iq,zq(ij,l,iq)
579c            print *,'q-->',ij,l,iq,q(ij,l,iq)
580             q(ij,l,iq)=zq(ij,l,iq)
581           ENDDO
582        ENDDO
583c$OMP END DO NOWAIT   
584      write(*,*) 'vlspltgen_loc 575'     
585
586c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
587        DO l=1,llm
588           DO ij=ijb,ije-iip1+1,iip1
589              q(ij+iim,l,iq)=q(ij,l,iq)
590           ENDDO
591        ENDDO
592c$OMP END DO NOWAIT 
593      write(*,*) 'vlspltgen_loc 583' 
594      ENDDO !DO iq=1,nqtot
595       
596      if (ok_iso_verif) then
597           call check_isotopes(q,ij_begin,ij_end,'vlspltgen_loc 557')
598      endif !if (ok_iso_verif) then
599
600c$OMP BARRIER
601
602cc$OMP MASTER     
603c      call WaitSendRequest(MyRequest1)
604c      call WaitSendRequest(MyRequest2)
605cc$OMP END MASTER
606cc$OMP BARRIER
607
608      write(*,*) 'vlspltgen 597: sortie' 
609      RETURN
610      END
Note: See TracBrowser for help on using the repository browser.