source: LMDZ6/trunk/libf/dyn3dmem/vlspltgen_loc.F @ 4500

Last change on this file since 4500 was 4469, checked in by Laurent Fairhead, 20 months ago

Replaced STOP instructions by calls to abort_gcm for a cleaner exit

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