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

Last change on this file since 2600 was 2600, checked in by Ehouarn Millour, 8 years ago

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