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

Last change on this file since 3093 was 2296, checked in by mvals, 5 years ago

Mars GCM:
dyn3dpar: Implementation in the parallel version of the dynamics of the dynamical transport of the isotopes using the same scheme as the one
implemented by Camille Risi in the Earth GCM (see LMDZ6/libf/dyn3dmem, and Risi et al. 2009: genealogy of the tracers+transport of the isotopic
ratio). This is added to prepare the future implementation of the HDO cycle in the GCM. These changes had been already added in the sequential part.
dyn3d: corrections to prevent from dividing by zero in the case of the use of the isotopes scheme.
traceur.def.isotopes: example of how to write the traceur.def if you want to use the new dynamics with the genealogy of the isotopes: air is the
father of all, H2O is the father of HDO.
MV
File size: 13.6 KB
Line 
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   --------------------------------------------------------------------
24      USE parallel_lmdz
25      USE mod_hallo
26      USE Write_Field_p
27      USE VAMPIR
28      USE infotrac, ONLY : nqtot,nqperes,nqdesc,nqfils,iqfils,
29     &    ok_iso_verif
30      USE comconst_mod, ONLY: cpp
31      IMPLICIT NONE
32
33c
34#include "dimensions.h"
35#include "paramet.h"
36
37c
38c   Arguments:
39c   ----------
40      INTEGER iadv(nqtot)
41      REAL masse(ip1jmp1,llm),pente_max
42      REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
43      REAL q(ip1jmp1,llm,nqtot)
44      REAL w(ip1jmp1,llm),pdt
45      REAL p(ip1jmp1,llmp1),teta(ip1jmp1,llm),pk(ip1jmp1,llm)
46c
47c      Local
48c   ---------
49c
50      INTEGER ij,l
51c
52      REAL,SAVE :: qsat(ip1jmp1,llm)
53      REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: zm
54      REAL,SAVE :: mu(ip1jmp1,llm)
55      REAL,SAVE :: mv(ip1jm,llm)
56      !REAL,SAVE :: mw(ip1jmp1,llm+1)
57      REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: mw !MVals
58      REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: zq
59      REAL zzpbar, zzw
60
61      REAL qmin,qmax
62      DATA qmin,qmax/0.,1.e33/
63
64c--pour rapport de melange saturant--
65
66      REAL rtt,retv,r2es,r3les,r3ies,r4les,r4ies,play
67      REAL ptarg,pdelarg,foeew,zdelta
68      REAL tempe(ip1jmp1)
69      INTEGER ijb,ije,iq,iq2,ifils
70      LOGICAL, SAVE :: firstcall=.TRUE.
71!$OMP THREADPRIVATE(firstcall)
72      type(request) :: MyRequest1
73      type(request) :: MyRequest2
74
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!$OMP MASTER
94            ALLOCATE(zm(ip1jmp1,llm,nqtot))
95            ALLOCATE(mw(ip1jmp1,llm+1,nqtot)) !MVals
96            ALLOCATE(zq(ip1jmp1,llm,nqtot))
97!$OMP END MASTER
98!$OMP BARRIER
99         END IF
100c-- Calcul de Qsat en chaque point
101c-- approximation: au milieu des couches play(l)=(p(l)+p(l+1))/2
102c   pour eviter une exponentielle.
103
104      call SetTag(MyRequest1,100)
105      call SetTag(MyRequest2,101)
106
107       
108        ijb=ij_begin-iip1
109        ije=ij_end+iip1
110        if (pole_nord) ijb=ij_begin
111        if (pole_sud) ije=ij_end
112       
113c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
114        DO l = 1, llm
115         DO ij = ijb, ije
116          tempe(ij) = teta(ij,l) * pk(ij,l) /cpp
117         ENDDO
118         DO ij = ijb, ije
119          zdelta = MAX( 0., SIGN(1., rtt - tempe(ij)) )
120          play   = 0.5*(p(ij,l)+p(ij,l+1))
121          qsat(ij,l) = MIN(0.5, r2es* FOEEW(tempe(ij),zdelta) / play )
122          qsat(ij,l) = qsat(ij,l) / ( 1. - retv * qsat(ij,l) )
123         ENDDO
124        ENDDO
125c$OMP END DO NOWAIT
126c      PRINT*,'Debut vlsplt version debug sans vlyqs'
127
128        zzpbar = 0.5 * pdt
129        zzw    = pdt
130
131      ijb=ij_begin
132      ije=ij_end
133      if (pole_nord) ijb=ijb+iip1
134      if (pole_sud)  ije=ije-iip1
135
136c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
137      DO l=1,llm
138        DO ij = ijb,ije
139            mu(ij,l)=pbaru(ij,l) * zzpbar
140         ENDDO
141      ENDDO
142c$OMP END DO NOWAIT
143     
144      ijb=ij_begin-iip1
145      ije=ij_end
146      if (pole_nord) ijb=ij_begin
147      if (pole_sud)  ije=ij_end-iip1
148
149c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
150      DO l=1,llm
151         DO ij=ijb,ije
152            mv(ij,l)=pbarv(ij,l) * zzpbar
153         ENDDO
154      ENDDO
155c$OMP END DO NOWAIT
156
157      ijb=ij_begin
158      ije=ij_end
159
160      DO iq=1,nqtot
161c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
162      DO l=1,llm
163         DO ij=ijb,ije
164            mw(ij,l,iq)=w(ij,l) * zzw
165         ENDDO
166      ENDDO
167c$OMP END DO NOWAIT
168      ENDDO
169
170      DO iq=1,nqtot
171c$OMP MASTER
172      DO ij=ijb,ije
173         mw(ij,llm+1,iq)=0.
174      ENDDO
175c$OMP END MASTER
176      ENDDO
177
178c      CALL SCOPY(ijp1llm,q,1,zq,1)
179c      CALL SCOPY(ijp1llm,masse,1,zm,1)
180
181       ijb=ij_begin
182       ije=ij_end
183
184      DO iq=1,nqtot
185c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
186        DO l=1,llm
187          zq(ijb:ije,l,iq)=q(ijb:ije,l,iq)
188          zm(ijb:ije,l,iq)=masse(ijb:ije,l)
189        ENDDO
190c$OMP END DO NOWAIT
191      ENDDO
192
193      ! verif temporaire
194      ijb=ij_begin
195      ije=ij_end
196      if (ok_iso_verif) then
197        call check_isotopes(zq,ijb,ije,'vlspltgen_loc 197')
198      endif !if (ok_iso_verif) then
199
200c$OMP BARRIER           
201!      DO iq=1,nqtot
202      DO iq=1,nqperes ! CRisi: on ne boucle que sur les pères= ceux qui sont transportés directement par l'air
203        if(iadv(iq) == 0) then
204       
205          cycle
206       
207        else if (iadv(iq)==10) then
208#ifdef _ADV_HALO     
209          call vlx_p(zq,pente_max,zm,mu,
210     &               ij_begin,ij_begin+2*iip1-1,iq)
211          call vlx_p(zq,pente_max,zm,mu,
212     &               ij_end-2*iip1+1,ij_end,iq)
213#else
214          call vlx_p(zq,pente_max,zm,mu,
215     &               ij_begin,ij_end,iq)
216#endif
217
218c$OMP MASTER
219          call VTb(VTHallo)
220c$OMP END MASTER
221          call Register_Hallo(zq(1,1,iq),ip1jmp1,llm,2,2,2,2,MyRequest1)
222          call Register_Hallo(zm(1,1,iq),ip1jmp1,llm,1,1,1,1,MyRequest1)
223! CRisi
224          do ifils=1,nqdesc(iq)
225            iq2=iqfils(ifils,iq)
226         call Register_Hallo(zq(1,1,iq2),ip1jmp1,llm,2,2,2,2,MyRequest1)
227         call Register_Hallo(zm(1,1,iq2),ip1jmp1,llm,1,1,1,1,MyRequest1)
228          enddo
229c$OMP MASTER
230          call VTe(VTHallo)
231c$OMP END MASTER
232        else if (iadv(iq)==14) then
233
234#ifdef _ADV_HALO           
235          call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat,
236     &                 ij_begin,ij_begin+2*iip1-1)
237          call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat,
238     &                 ij_end-2*iip1+1,ij_end)
239#else
240
241          call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat,
242     &                 ij_begin,ij_end)
243#endif
244
245c$OMP MASTER
246          call VTb(VTHallo)
247c$OMP END MASTER
248
249          call Register_Hallo(zq(1,1,iq),ip1jmp1,llm,2,2,2,2,MyRequest1)
250          call Register_Hallo(zm(1,1,iq),ip1jmp1,llm,1,1,1,1,MyRequest1)
251!CRisi
252          do ifils=1,nqdesc(iq)
253            iq2=iqfils(ifils,iq)
254         call Register_Hallo(zq(1,1,iq2),ip1jmp1,llm,2,2,2,2,MyRequest1)
255         call Register_Hallo(zm(1,1,iq2),ip1jmp1,llm,1,1,1,1,MyRequest1)
256          enddo
257c$OMP MASTER
258          call VTe(VTHallo)
259c$OMP END MASTER
260        else
261       
262          stop 'vlspltgen_p : schema non parallelise'
263     
264        endif
265     
266      enddo !DO iq=1,nqperes
267     
268     
269c$OMP BARRIER     
270c$OMP MASTER     
271      call VTb(VTHallo)
272c$OMP END MASTER
273
274      call SendRequest(MyRequest1)
275
276c$OMP MASTER
277      call VTe(VTHallo)
278c$OMP END MASTER       
279c$OMP BARRIER
280
281!      ! verif temporaire
282!      ijb=ij_begin-2*iip1
283!      ije=ij_end+2*iip1 
284!      if (pole_nord) ijb=ij_begin
285!      if (pole_sud)  ije=ij_end 
286      if (ok_iso_verif) then
287           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 293')
288      endif !if (ok_iso_verif) then
289
290!      do iq=1,nqtot
291      do iq=1,nqperes !CRisi
292
293        if(iadv(iq) == 0) then
294       
295          cycle
296       
297        else if (iadv(iq)==10) then
298
299#ifdef _ADV_HALLO
300          call vlx_p(zq,pente_max,zm,mu,
301     &               ij_begin+2*iip1,ij_end-2*iip1,iq)
302#endif       
303        else if (iadv(iq)==14) then
304#ifdef _ADV_HALLO
305          call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat,
306     &                 ij_begin+2*iip1,ij_end-2*iip1)
307#endif   
308        else
309       
310          stop 'vlspltgen_p : schema non parallelise'
311     
312        endif
313     
314      enddo
315c$OMP BARRIER     
316c$OMP MASTER
317      call VTb(VTHallo)
318c$OMP END MASTER
319
320!      call WaitRecvRequest(MyRequest1)
321!      call WaitSendRequest(MyRequest1)
322c$OMP BARRIER
323       call WaitRequest(MyRequest1)
324
325
326c$OMP MASTER
327      call VTe(VTHallo)
328c$OMP END MASTER
329c$OMP BARRIER
330
331      if (ok_iso_verif) then
332           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 326')
333      endif !if (ok_iso_verif) then       
334      if (ok_iso_verif) then
335           ijb=ij_begin-2*iip1
336           ije=ij_end+2*iip1
337           if (pole_nord) ijb=ij_begin
338           if (pole_sud)  ije=ij_end
339           call check_isotopes(zq,ijb,ije,'vlspltgen_loc 336')
340      endif !if (ok_iso_verif) then
341
342!      do iq=1,nqtot
343      do iq=1,nqperes !CRisi
344
345        if(iadv(iq) == 0) then
346       
347          cycle
348       
349        else if (iadv(iq)==10) then
350       
351          call vly_p(zq,pente_max,zm,mv,iq)
352 
353        else if (iadv(iq)==14) then
354     
355          call vlyqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv,qsat)
356 
357        else
358       
359          stop 'vlspltgen_p : schema non parallelise'
360     
361        endif
362       
363       enddo
364
365      if (ok_iso_verif) then
366           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 357')
367      endif !if (ok_iso_verif) then
368
369!      do iq=1,nqtot
370      do iq=1,nqperes !CRisi
371
372        if(iadv(iq) == 0) then
373         
374          cycle
375       
376        else if (iadv(iq)==10 .or. iadv(iq)==14 ) then
377
378c$OMP BARRIER       
379#ifdef _ADV_HALLO
380          call vlz_p(zq,pente_max,zm,mw,
381     &               ij_begin,ij_begin+2*iip1-1,iq)
382          call vlz_p(zq,pente_max,zm,mw,
383     &               ij_end-2*iip1+1,ij_end,iq)
384#else
385          call vlz_p(zq,pente_max,zm,mw,
386     &               ij_begin,ij_end,iq)
387#endif
388c$OMP BARRIER
389
390c$OMP MASTER
391          call VTb(VTHallo)
392c$OMP END MASTER
393
394          call Register_Hallo(zq(1,1,iq),ip1jmp1,llm,2,2,2,2,MyRequest2)
395          call Register_Hallo(zm(1,1,iq),ip1jmp1,llm,1,1,1,1,MyRequest2)
396          ! CRisi
397          do ifils=1,nqdesc(iq)
398            iq2=iqfils(ifils,iq)
399         call Register_Hallo(zq(1,1,iq2),ip1jmp1,llm,2,2,2,2,MyRequest2)
400         call Register_Hallo(zm(1,1,iq2),ip1jmp1,llm,1,1,1,1,MyRequest2)
401          enddo
402c$OMP MASTER
403          call VTe(VTHallo)
404c$OMP END MASTER       
405c$OMP BARRIER
406        else
407       
408          stop 'vlspltgen_p : schema non parallelise'
409     
410        endif
411     
412      enddo
413c$OMP BARRIER     
414
415c$OMP MASTER       
416      call VTb(VTHallo)
417c$OMP END MASTER
418
419      call SendRequest(MyRequest2)
420
421c$OMP MASTER
422      call VTe(VTHallo)
423c$OMP END MASTER       
424
425      if (ok_iso_verif) then
426           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 429')
427      endif !if (ok_iso_verif) then
428
429c$OMP BARRIER
430!      do iq=1,nqtot
431      do iq=1,nqperes !CRisi
432
433        if(iadv(iq) == 0) then
434         
435          cycle
436       
437        else if (iadv(iq)==10 .or. iadv(iq)==14 ) then
438c$OMP BARRIER       
439
440#ifdef _ADV_HALLO
441          call vlz_p(zq,pente_max,zm,mw,
442     &               ij_begin+2*iip1,ij_end-2*iip1,iq)
443#endif
444
445c$OMP BARRIER       
446        else
447       
448          stop 'vlspltgen_p : schema non parallelise'
449     
450        endif
451     
452      enddo
453
454c$OMP BARRIER
455c$OMP MASTER
456      call VTb(VTHallo)
457c$OMP END MASTER
458
459!      call WaitRecvRequest(MyRequest2)
460!      call WaitSendRequest(MyRequest2)
461c$OMP BARRIER
462       CALL WaitRequest(MyRequest2)
463
464c$OMP MASTER
465      call VTe(VTHallo)
466c$OMP END MASTER
467c$OMP BARRIER
468
469      !write(*,*) 'vlspltgen_loc 494'
470      if (ok_iso_verif) then
471           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 461')
472      endif !if (ok_iso_verif) then
473
474!      do iq=1,nqtot
475      do iq=1,nqperes !CRisi
476
477        if(iadv(iq) == 0) then
478       
479          cycle
480       
481        else if (iadv(iq)==10) then
482       
483          call vly_p(zq,pente_max,zm,mv,iq)
484 
485        else if (iadv(iq)==14) then
486     
487          call vlyqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv,qsat)
488 
489        else
490       
491          stop 'vlspltgen_p : schema non parallelise'
492     
493        endif
494       
495       enddo !do iq=1,nqperes
496
497      if (ok_iso_verif) then
498           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 493')
499      endif !if (ok_iso_verif) then
500
501!      do iq=1,nqtot
502      do iq=1,nqperes !CRisi
503
504        if(iadv(iq) == 0) then
505         
506          cycle
507       
508        else if (iadv(iq)==10) then
509       
510          call vlx_p(zq,pente_max,zm,mu,
511     &               ij_begin,ij_end,iq)
512 
513        else if (iadv(iq)==14) then
514     
515          call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat,
516     &                 ij_begin,ij_end)
517 
518        else
519       
520          stop 'vlspltgen_p : schema non parallelise'
521     
522        endif
523       
524       enddo !do iq=1,nqperes
525
526      !write(*,*) 'vlspltgen 550: apres derniere serie de call vlx'
527      if (ok_iso_verif) then
528           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 521')
529      endif !if (ok_iso_verif) then
530     
531      ijb=ij_begin
532      ije=ij_end
533c$OMP BARRIER     
534
535
536      DO iq=1,nqtot
537
538c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
539        DO l=1,llm
540           DO ij=ijb,ije
541c             print *,'zq-->',ij,l,iq,zq(ij,l,iq)
542c            print *,'q-->',ij,l,iq,q(ij,l,iq)
543             q(ij,l,iq)=zq(ij,l,iq)
544           ENDDO
545        ENDDO
546c$OMP END DO NOWAIT         
547
548c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
549        DO l=1,llm
550           DO ij=ijb,ije-iip1+1,iip1
551              q(ij+iim,l,iq)=q(ij,l,iq)
552           ENDDO
553        ENDDO
554c$OMP END DO NOWAIT 
555
556      ENDDO
557       
558      if (ok_iso_verif) then
559           call check_isotopes(q,ij_begin,ij_end,'vlspltgen_loc 557')
560      endif !if (ok_iso_verif) then
561
562c$OMP BARRIER
563
564cc$OMP MASTER     
565c      call WaitSendRequest(MyRequest1)
566c      call WaitSendRequest(MyRequest2)
567cc$OMP END MASTER
568cc$OMP BARRIER
569
570      RETURN
571      END
Note: See TracBrowser for help on using the repository browser.