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

Last change on this file since 3990 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
File size: 15.2 KB
RevLine 
[1632]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   --------------------------------------------------------------------
[1864]25      USE parallel_lmdz
[1632]26      USE mod_hallo
27      USE Write_Field_loc
28      USE VAMPIR
[2298]29      ! CRisi: on rajoute variables utiles d'infotrac 
30      USE infotrac, ONLY : nqtot,nqperes,nqdesc,nqfils,iqfils,
31     &    ok_iso_verif
[1632]32      USE vlspltgen_mod
[2641]33      USE comconst_mod, ONLY: cpp
[1632]34      IMPLICIT NONE
35
36c
[2641]37      include "dimensions.h"
38      include "paramet.h"
[1632]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)
[2298]66      INTEGER ijb,ije,iq,iq2,ifils
[1632]67      LOGICAL, SAVE :: firstcall=.TRUE.
68!$OMP THREADPRIVATE(firstcall)
[1864]69      type(request),SAVE :: MyRequest1
70!$OMP THREADPRIVATE(MyRequest1)
71      type(request),SAVE :: MyRequest2
72!$OMP THREADPRIVATE(MyRequest2)
[1632]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       
[2641]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       
[1632]105c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[2641]106        DO l = 1, llm
[1632]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
[2298]152      DO iq=1,nqtot
[1632]153c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
154      DO l=1,llm
155         DO ij=ijb,ije
[2298]156            mw(ij,l,iq)=w(ij,l) * zzw
[1632]157         ENDDO
158      ENDDO
159c$OMP END DO NOWAIT
[2298]160      ENDDO
[1632]161
[2298]162      DO iq=1,nqtot 
[1632]163c$OMP MASTER
164      DO ij=ijb,ije
[2298]165         mw(ij,llm+1,iq)=0.
[1632]166      ENDDO
167c$OMP END MASTER
[2298]168      ENDDO
[1632]169
170c      CALL SCOPY(ijp1llm,q,1,zq,1)
171c      CALL SCOPY(ijp1llm,masse,1,zm,1)
172
173       ijb=ij_begin
174       ije=ij_end
175
[2298]176      DO iq=1,nqtot       
[1632]177c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
178        DO l=1,llm
179          zq(ijb:ije,l,iq)=q(ijb:ije,l,iq)
180          zm(ijb:ije,l,iq)=masse(ijb:ije,l)
181        ENDDO
182c$OMP END DO NOWAIT
183      ENDDO
184
[2298]185#ifdef DEBUG_IO     
[1632]186       CALL WriteField_u('mu',mu)
187       CALL WriteField_v('mv',mv)
188       CALL WriteField_u('mw',mw)
189       CALL WriteField_u('qsat',qsat)
190#endif
191
[2298]192      ! verif temporaire
193      ijb=ij_begin
194      ije=ij_end 
195      if (ok_iso_verif) then
196        call check_isotopes(zq,ijb,ije,'vlspltgen_loc 191')
197      endif !if (ok_iso_verif) then   
198
[1632]199c$OMP BARRIER           
[2298]200!      DO iq=1,nqtot
201      DO iq=1,nqperes ! CRisi: on ne boucle que sur les pères= ceux qui sont transportés directement par l'air
202       !write(*,*) 'vlspltgen 192: iq,iadv=',iq,iadv(iq)
[1632]203#ifdef DEBUG_IO   
204       CALL WriteField_u('zq',zq(:,:,iq))
[2298]205       CALL WriteField_u('zm',zm(:,:,iq))
[1632]206#endif
207        if(iadv(iq) == 0) then
[2641]208       
209          cycle
210       
211        else if (iadv(iq)==10) then
[1632]212
213#ifdef _ADV_HALO       
[2298]214! CRisi: on ajoute les nombres de fils et tableaux des fils
215! On suppose qu'on ne peut advecter les fils que par le schéma 10. 
216          call vlx_loc(zq,pente_max,zm,mu,
[2641]217     &                     ij_begin,ij_begin+2*iip1-1,iq)
218          call vlx_loc(zq,pente_max,zm,mu,
[2298]219     &               ij_end-2*iip1+1,ij_end,iq)
[1632]220#else
[2641]221          call vlx_loc(zq,pente_max,zm,mu,
222     &                     ij_begin,ij_end,iq)
[1632]223#endif
224
225c$OMP MASTER
226          call VTb(VTHallo)
227c$OMP END MASTER
228          call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest1)
229          call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest1)
[2298]230! CRisi
231          do ifils=1,nqdesc(iq)
232            iq2=iqfils(ifils,iq)
233            call Register_Hallo_u(zq(:,:,iq2),llm,2,2,2,2,MyRequest1)
234            call Register_Hallo_u(zm(:,:,iq2),llm,1,1,1,1,MyRequest1)
235          enddo
[1632]236
237c$OMP MASTER
238          call VTe(VTHallo)
239c$OMP END MASTER
[2641]240        else if (iadv(iq)==14) then
[1632]241
242#ifdef _ADV_HALO           
[2298]243          call vlxqs_loc(zq,pente_max,zm,mu,
244     &                   qsat,ij_begin,ij_begin+2*iip1-1,iq)
245          call vlxqs_loc(zq,pente_max,zm,mu,
246     &                   qsat,ij_end-2*iip1+1,ij_end,iq)
[1632]247#else
[2298]248          call vlxqs_loc(zq,pente_max,zm,mu,
249     &                   qsat,ij_begin,ij_end,iq)
[1632]250#endif
251
252c$OMP MASTER
253          call VTb(VTHallo)
254c$OMP END MASTER
255
256          call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest1)
257          call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest1)
[2298]258          do ifils=1,nqdesc(iq)
259            iq2=iqfils(ifils,iq)
260            call Register_Hallo_u(zq(:,:,iq2),llm,2,2,2,2,MyRequest1)
261            call Register_Hallo_u(zm(:,:,iq2),llm,1,1,1,1,MyRequest1)
262          enddo
[1632]263
264c$OMP MASTER
265          call VTe(VTHallo)
266c$OMP END MASTER
267        else
[2641]268       
269          stop 'vlspltgen_p : schema non parallelise'
[1632]270     
271        endif
272     
[2298]273      enddo !DO iq=1,nqperes
[1632]274     
275     
276c$OMP BARRIER     
277c$OMP MASTER     
278      call VTb(VTHallo)
279c$OMP END MASTER
280
281      call SendRequest(MyRequest1)
282
283c$OMP MASTER
284      call VTe(VTHallo)
285c$OMP END MASTER       
286c$OMP BARRIER
287
[2298]288      ! verif temporaire
289      ijb=ij_begin-2*iip1
290      ije=ij_end+2*iip1 
291      if (pole_nord) ijb=ij_begin
292      if (pole_sud)  ije=ij_end 
293      if (ok_iso_verif) then
294           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 280')
295      endif !if (ok_iso_verif) then
296
297      do iq=1,nqperes
298        !write(*,*) 'vlspltgen 279: iq=',iq
299
[1632]300        if(iadv(iq) == 0) then
[2641]301       
302          cycle
303       
304        else if (iadv(iq)==10) then
[1632]305
306#ifdef _ADV_HALLO
[2298]307          call vlx_loc(zq,pente_max,zm,mu,
308     &                 ij_begin+2*iip1,ij_end-2*iip1,iq)
[1632]309#endif       
[2641]310        else if (iadv(iq)==14) then
[1632]311#ifdef _ADV_HALLO
[2298]312          call vlxqs_loc(zq,pente_max,zm,mu,
313     &                    qsat,ij_begin+2*iip1,ij_end-2*iip1,iq)
[1632]314#endif   
315        else
[2641]316       
317          stop 'vlspltgen_p : schema non parallelise'
[1632]318     
319        endif
320     
321      enddo
322c$OMP BARRIER     
323c$OMP MASTER
324      call VTb(VTHallo)
325c$OMP END MASTER
326
327!      call WaitRecvRequest(MyRequest1)
328!      call WaitSendRequest(MyRequest1)
329c$OMP BARRIER
330       call WaitRequest(MyRequest1)
331
332
333c$OMP MASTER
334      call VTe(VTHallo)
335c$OMP END MASTER
336c$OMP BARRIER
337
[2298]338     
339      if (ok_iso_verif) then
340           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 326')
341      endif !if (ok_iso_verif) then       
342      if (ok_iso_verif) then
343           ijb=ij_begin-2*iip1
344           ije=ij_end+2*iip1
345           if (pole_nord) ijb=ij_begin
346           if (pole_sud)  ije=ij_end
347           call check_isotopes(zq,ijb,ije,'vlspltgen_loc 336')
348      endif !if (ok_iso_verif) then 
[1632]349
[2298]350      do iq=1,nqperes
351       !write(*,*) 'vlspltgen 321: iq=',iq
[1632]352#ifdef DEBUG_IO   
353       CALL WriteField_u('zq',zq(:,:,iq))
354       CALL WriteField_u('zm',zm(:,:,iq))
355#endif
[2298]356
[1632]357        if(iadv(iq) == 0) then
358       
[2641]359          cycle
360       
361        else if (iadv(iq)==10) then
362       
[2298]363          call vly_loc(zq,pente_max,zm,mv,iq)
[1632]364 
[2641]365        else if (iadv(iq)==14) then
[1632]366     
[2298]367          call vlyqs_loc(zq,pente_max,zm,mv,
368     &                   qsat,iq)
[1632]369 
370        else
[2641]371       
372          stop 'vlspltgen_p : schema non parallelise'
[1632]373     
374        endif
375       
376       enddo
377
[2298]378      if (ok_iso_verif) then
379           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 357')
380      endif !if (ok_iso_verif) then
[1632]381
[2298]382      do iq=1,nqperes
383      !write(*,*) 'vlspltgen 349: iq=',iq
[1632]384#ifdef DEBUG_IO   
385       CALL WriteField_u('zq',zq(:,:,iq))
386       CALL WriteField_u('zm',zm(:,:,iq))
387#endif
388        if(iadv(iq) == 0) then
[2641]389         
390          cycle
391       
392        else if (iadv(iq)==10 .or. iadv(iq)==14 ) then
[1632]393
394c$OMP BARRIER       
395#ifdef _ADV_HALLO
[2298]396          call vlz_loc(zq,pente_max,zm,mw,
397     &               ij_begin,ij_begin+2*iip1-1,iq)
398          call vlz_loc(zq,pente_max,zm,mw,
399     &               ij_end-2*iip1+1,ij_end,iq)
[1632]400#else
[2298]401          call vlz_loc(zq,pente_max,zm,mw,
402     &               ij_begin,ij_end,iq)
[1632]403#endif
404c$OMP BARRIER
405
406c$OMP MASTER
407          call VTb(VTHallo)
408c$OMP END MASTER
409
410          call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest2)
411          call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest2)
[2298]412          ! CRisi
413          do ifils=1,nqdesc(iq)
414            iq2=iqfils(ifils,iq)
415            call Register_Hallo_u(zq(:,:,iq2),llm,2,2,2,2,MyRequest2)
416            call Register_Hallo_u(zm(:,:,iq2),llm,1,1,1,1,MyRequest2)
417          enddo     
[1632]418c$OMP MASTER
419          call VTe(VTHallo)
[2641]420c$OMP END MASTER       
[1632]421c$OMP BARRIER
422        else
[2641]423       
424          stop 'vlspltgen_p : schema non parallelise'
[1632]425     
426        endif
427     
428      enddo
429c$OMP BARRIER     
430
431c$OMP MASTER       
432      call VTb(VTHallo)
433c$OMP END MASTER
434
435      call SendRequest(MyRequest2)
436
437c$OMP MASTER
438      call VTe(VTHallo)
[2641]439c$OMP END MASTER       
[1632]440
[2298]441
442      if (ok_iso_verif) then
443           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 429')
444      endif !if (ok_iso_verif) then
445
[1632]446c$OMP BARRIER
[2298]447      do iq=1,nqperes
448      !write(*,*) 'vlspltgen 409: iq=',iq
[1632]449
450        if(iadv(iq) == 0) then
[2641]451         
452          cycle
453       
454        else if (iadv(iq)==10 .or. iadv(iq)==14 ) then
[1632]455c$OMP BARRIER       
456
457#ifdef _ADV_HALLO
[2298]458          call vlz_loc(zq,pente_max,zm,mw,
459     &               ij_begin+2*iip1,ij_end-2*iip1,iq)
[1632]460#endif
461
462c$OMP BARRIER       
463        else
[2641]464       
465          stop 'vlspltgen_p : schema non parallelise'
[1632]466     
467        endif
468     
469      enddo
[2298]470      !write(*,*) 'vlspltgen_loc 476'
[1632]471
472c$OMP BARRIER
[2298]473      !write(*,*) 'vlspltgen_loc 477'
[1632]474c$OMP MASTER
475      call VTb(VTHallo)
476c$OMP END MASTER
477
478!      call WaitRecvRequest(MyRequest2)
479!      call WaitSendRequest(MyRequest2)
480c$OMP BARRIER
481       CALL WaitRequest(MyRequest2)
482
483c$OMP MASTER
484      call VTe(VTHallo)
485c$OMP END MASTER
486c$OMP BARRIER
487
488
[2298]489      !write(*,*) 'vlspltgen_loc 494'
490      if (ok_iso_verif) then
491           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 461')
492      endif !if (ok_iso_verif) then
493
494      do iq=1,nqperes
495      !write(*,*) 'vlspltgen 449: iq=',iq
[1632]496#ifdef DEBUG_IO   
497       CALL WriteField_u('zq',zq(:,:,iq))
498       CALL WriteField_u('zm',zm(:,:,iq))
499#endif
500        if(iadv(iq) == 0) then
501       
[2641]502          cycle
503       
504        else if (iadv(iq)==10) then
505       
[2298]506          call vly_loc(zq,pente_max,zm,mv,iq)
[1632]507 
[2641]508        else if (iadv(iq)==14) then
[1632]509     
[2298]510          call vlyqs_loc(zq,pente_max,zm,mv,
511     &                   qsat,iq)
[1632]512 
513        else
[2641]514       
515          stop 'vlspltgen_p : schema non parallelise'
[1632]516     
517        endif
518       
[2298]519       enddo !do iq=1,nqperes
[1632]520
[2298]521      if (ok_iso_verif) then
522           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 493')
523      endif !if (ok_iso_verif) then
[1632]524
[2298]525      do iq=1,nqperes
526      !write(*,*) 'vlspltgen 477: iq=',iq
[1632]527#ifdef DEBUG_IO   
528       CALL WriteField_u('zq',zq(:,:,iq))
529       CALL WriteField_u('zm',zm(:,:,iq))
530#endif
531        if(iadv(iq) == 0) then
[2641]532         
533          cycle
[1632]534       
[2641]535        else if (iadv(iq)==10) then
536       
[2298]537          call vlx_loc(zq,pente_max,zm,mu,
538     &               ij_begin,ij_end,iq)
[1632]539 
[2641]540        else if (iadv(iq)==14) then
[1632]541     
[2298]542          call vlxqs_loc(zq,pente_max,zm,mu,
543     &                 qsat, ij_begin,ij_end,iq)
[1632]544 
545        else
[2641]546       
[1632]547          stop 'vlspltgen_p : schema non parallelise'
548     
549        endif
550       
[2298]551       enddo !do iq=1,nqperes
[1632]552
[2298]553      !write(*,*) 'vlspltgen 550: apres derniere serie de call vlx'
554      if (ok_iso_verif) then
555           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 521')
556      endif !if (ok_iso_verif) then
[1632]557     
558      ijb=ij_begin
559      ije=ij_end
[2298]560      !write(*,*) 'vlspltgen_loc 557'
[1632]561c$OMP BARRIER     
562
[2298]563      !write(*,*) 'vlspltgen_loc 559' 
[1632]564      DO iq=1,nqtot
[2298]565       !write(*,*) 'vlspltgen_loc 561, iq=',iq 
[1632]566#ifdef DEBUG_IO   
567       CALL WriteField_u('zq',zq(:,:,iq))
568       CALL WriteField_u('zm',zm(:,:,iq))
569#endif
570c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
571        DO l=1,llm
572           DO ij=ijb,ije
573c             print *,'zq-->',ij,l,iq,zq(ij,l,iq)
[2641]574c             print *,'q-->',ij,l,iq,q(ij,l,iq)
575             q(ij,l,iq)=zq(ij,l,iq)
[1632]576           ENDDO
577        ENDDO
[2298]578c$OMP END DO NOWAIT   
579      !write(*,*) 'vlspltgen_loc 575'     
[1632]580
581c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
582        DO l=1,llm
583           DO ij=ijb,ije-iip1+1,iip1
584              q(ij+iim,l,iq)=q(ij,l,iq)
585           ENDDO
586        ENDDO
587c$OMP END DO NOWAIT 
[2298]588      !write(*,*) 'vlspltgen_loc 583' 
589      ENDDO !DO iq=1,nqtot
590       
591      if (ok_iso_verif) then
592           call check_isotopes(q,ij_begin,ij_end,'vlspltgen_loc 557')
593      endif !if (ok_iso_verif) then
[1632]594
595c$OMP BARRIER
596
597cc$OMP MASTER     
598c      call WaitSendRequest(MyRequest1)
599c      call WaitSendRequest(MyRequest2)
600cc$OMP END MASTER
601cc$OMP BARRIER
602
[2298]603      !write(*,*) 'vlspltgen 597: sortie' 
[1632]604      RETURN
605      END
Note: See TracBrowser for help on using the repository browser.