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

Last change on this file since 2361 was 2286, checked in by Ehouarn Millour, 10 years ago
  • fix in dyn3d the array out of bound issue (in qminimum) that was only corrected in dyn3dmem in rev 2285.
  • comment out many invasive debug writes.

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