source: LMDZ6/branches/LMDZ_ECRad/libf/dyn3dmem/vlspltgen_loc.F @ 5731

Last change on this file since 5731 was 4482, checked in by lguez, 3 years ago

Sync latest trunk changes to branch LMDZ_ECRad
  • 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
Line 
1
2!     
3! $Header$
4!
5       SUBROUTINE vlspltgen_loc( q,pente_max,masse,w,pbaru,pbarv,
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   --------------------------------------------------------------------
26      USE parallel_lmdz
27      USE mod_hallo
28      USE Write_Field_loc
29      USE VAMPIR
30      ! CRisi: on rajoute variables utiles d'infotrac 
31      USE infotrac, ONLY : nqtot, tracers, isoCheck
32      USE vlspltgen_mod
33      USE comconst_mod, ONLY: cpp
34      IMPLICIT NONE
35
36c
37      include "dimensions.h"
38      include "paramet.h"
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)
65      INTEGER ijb,ije,iq,iq2,ifils
66      LOGICAL, SAVE :: firstcall=.TRUE.
67!$OMP THREADPRIVATE(firstcall)
68      type(request),SAVE :: MyRequest1
69!$OMP THREADPRIVATE(MyRequest1)
70      type(request),SAVE :: MyRequest2
71!$OMP THREADPRIVATE(MyRequest2)
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       
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       
104c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
105        DO l = 1, llm
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
151      DO iq=1,nqtot
152c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
153      DO l=1,llm
154         DO ij=ijb,ije
155            mw(ij,l,iq)=w(ij,l) * zzw
156         ENDDO
157      ENDDO
158c$OMP END DO NOWAIT
159      ENDDO
160
161      DO iq=1,nqtot 
162c$OMP MASTER
163      DO ij=ijb,ije
164         mw(ij,llm+1,iq)=0.
165      ENDDO
166c$OMP END MASTER
167      ENDDO
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
175      DO iq=1,nqtot       
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
184#ifdef DEBUG_IO     
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
191      ! verif temporaire
192      ijb=ij_begin
193      ije=ij_end 
194      call check_isotopes(zq,ijb,ije,'vlspltgen_loc 191')
195
196c$OMP BARRIER           
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
201#ifdef DEBUG_IO   
202        CALL WriteField_u('zq',zq(:,:,iq))
203        CALL WriteField_u('zm',zm(:,:,iq))
204#endif
205        SELECT CASE(tracers(iq)%iadv)
206          CASE(0); CYCLE
207          CASE(10)
208#ifdef _ADV_HALO       
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,
212     &                     ij_begin,ij_begin+2*iip1-1,iq)
213          call vlx_loc(zq,pente_max,zm,mu,
214     &               ij_end-2*iip1+1,ij_end,iq)
215#else
216          call vlx_loc(zq,pente_max,zm,mu,
217     &                     ij_begin,ij_end,iq)
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)
225! CRisi
226          do ifils=1,tracers(iq)%nqDescen
227            iq2=tracers(iq)%iqDescen(ifils)
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
231
232c$OMP MASTER
233          call VTe(VTHallo)
234c$OMP END MASTER
235          CASE(14)
236#ifdef _ADV_HALO           
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)
241#else
242          call vlxqs_loc(zq,pente_max,zm,mu,
243     &                   qsat,ij_begin,ij_end,iq)
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)
252          do ifils=1,tracers(iq)%nqDescen
253            iq2=tracers(iq)%iqDescen(ifils)
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
257
258c$OMP MASTER
259          call VTe(VTHallo)
260c$OMP END MASTER
261          CASE DEFAULT
262             CALL abort_gcm("vlspltgen_loc","schema non parallelise",1)
263        END SELECT
264     
265      enddo !DO iq=1,nqtot
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
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 
285      call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 280')
286
287      do iq=1,nqtot
288        IF(tracers(iq)%parent /= 'air') CYCLE
289        !write(*,*) 'vlspltgen 279: iq=',iq
290
291        SELECT CASE(tracers(iq)%iadv)
292          CASE(0); CYCLE
293          CASE(10)
294#ifdef _ADV_HALLO
295          call vlx_loc(zq,pente_max,zm,mu,
296     &                 ij_begin+2*iip1,ij_end-2*iip1,iq)
297#endif       
298          CASE(14)
299#ifdef _ADV_HALLO
300          call vlxqs_loc(zq,pente_max,zm,mu,
301     &                    qsat,ij_begin+2*iip1,ij_end-2*iip1,iq)
302#endif   
303          CASE DEFAULT
304          CALL abort_gcm("vlspltgen_p","schema non parallelise",1)
305        END SELECT
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
324     
325      IF(isoCheck) THEN
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')
332      END IF
333
334      do iq = 1, nqtot
335       IF(tracers(iq)%parent /= 'air') CYCLE
336       !write(*,*) 'vlspltgen 321: iq=',iq
337#ifdef DEBUG_IO   
338       CALL WriteField_u('zq',zq(:,:,iq))
339       CALL WriteField_u('zm',zm(:,:,iq))
340#endif
341
342        SELECT CASE(tracers(iq)%iadv)
343          CASE(0); CYCLE
344          CASE(10); call   vly_loc(zq,pente_max,zm,mv,     iq)
345          CASE(14); call vlyqs_loc(zq,pente_max,zm,mv,qsat,iq)
346          CASE DEFAULT
347          CALL abort_gcm("vlspltgen_p","schema non parallelise",1)
348        END SELECT
349       
350       enddo
351
352      call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 357')
353
354      do iq = 1, nqtot
355       IF(tracers(iq)%parent /= 'air') CYCLE
356      !write(*,*) 'vlspltgen 349: iq=',iq
357#ifdef DEBUG_IO   
358       CALL WriteField_u('zq',zq(:,:,iq))
359       CALL WriteField_u('zm',zm(:,:,iq))
360#endif
361        SELECT CASE(tracers(iq)%iadv)
362          CASE(0); CYCLE
363          CASE(10,14)
364c$OMP BARRIER       
365#ifdef _ADV_HALLO
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)
370#else
371          call vlz_loc(zq,pente_max,zm,mw,
372     &               ij_begin,ij_end,iq)
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)
382          ! CRisi
383          do ifils=1,tracers(iq)%nqDescen
384            iq2=tracers(iq)%iqDescen(ifils)
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     
388c$OMP MASTER
389          call VTe(VTHallo)
390c$OMP END MASTER       
391c$OMP BARRIER
392          CASE DEFAULT
393           
394            CALL abort_gcm("vlspltgen_p","schema non parallelise",1)
395         END SELECT
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)
408c$OMP END MASTER       
409
410
411      call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 429')
412
413c$OMP BARRIER
414      do iq=1,nqtot
415        IF(tracers(iq)%parent /= 'air') CYCLE
416      !write(*,*) 'vlspltgen 409: iq=',iq
417
418        SELECT CASE(tracers(iq)%iadv)
419          CASE(0); CYCLE
420          CASE(10,14)
421c$OMP BARRIER       
422
423#ifdef _ADV_HALLO
424          call vlz_loc(zq,pente_max,zm,mw,
425     &               ij_begin+2*iip1,ij_end-2*iip1,iq)
426#endif
427
428c$OMP BARRIER       
429          CASE DEFAULT
430          CALL abort_gcm("vlspltgen_p","schema non parallelise",1)
431        END SELECT
432     
433      enddo
434      !write(*,*) 'vlspltgen_loc 476'
435
436c$OMP BARRIER
437      !write(*,*) 'vlspltgen_loc 477'
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
453      !write(*,*) 'vlspltgen_loc 494'
454      call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 461')
455
456      do iq=1,nqtot
457        IF(tracers(iq)%parent /= 'air') CYCLE
458      !write(*,*) 'vlspltgen 449: iq=',iq
459#ifdef DEBUG_IO   
460       CALL WriteField_u('zq',zq(:,:,iq))
461       CALL WriteField_u('zm',zm(:,:,iq))
462#endif
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)
467          CASE DEFAULT
468             CALL abort_gcm("vlspltgen_p","schema non parallelise",1)
469        END SELECT
470       
471       enddo !do iq=1,nqtot
472
473      call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 493')
474
475      do iq=1,nqtot
476        IF(tracers(iq)%parent /= 'air') CYCLE
477      !write(*,*) 'vlspltgen 477: iq=',iq
478#ifdef DEBUG_IO   
479       CALL WriteField_u('zq',zq(:,:,iq))
480       CALL WriteField_u('zm',zm(:,:,iq))
481#endif
482        SELECT CASE(tracers(iq)%iadv)
483          CASE(0); CYCLE
484          CASE(10); call   vlx_loc(zq,pente_max,zm,mu,
485     &               ij_begin,ij_end,iq)
486          CASE(14); call vlxqs_loc(zq,pente_max,zm,mu,
487     &                 qsat, ij_begin,ij_end,iq)
488          CASE DEFAULT
489          CALL abort_gcm("vlspltgen_p","schema non parallelise",1)
490        END SELECT
491       
492       enddo !do iq=1,nqtot
493
494      !write(*,*) 'vlspltgen 550: apres derniere serie de call vlx'
495      call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 521')
496     
497      ijb=ij_begin
498      ije=ij_end
499      !write(*,*) 'vlspltgen_loc 557'
500c$OMP BARRIER     
501
502      !write(*,*) 'vlspltgen_loc 559' 
503      DO iq=1,nqtot
504       !write(*,*) 'vlspltgen_loc 561, iq=',iq 
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)
513c             print *,'q-->',ij,l,iq,q(ij,l,iq)
514             q(ij,l,iq)=zq(ij,l,iq)
515           ENDDO
516        ENDDO
517c$OMP END DO NOWAIT   
518      !write(*,*) 'vlspltgen_loc 575'     
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 
527      !write(*,*) 'vlspltgen_loc 583' 
528      ENDDO !DO iq=1,nqtot
529       
530      call check_isotopes(q,ij_begin,ij_end,'vlspltgen_loc 557')
531
532c$OMP BARRIER
533
534cc$OMP MASTER     
535c      call WaitSendRequest(MyRequest1)
536c      call WaitSendRequest(MyRequest2)
537cc$OMP END MASTER
538cc$OMP BARRIER
539
540      !write(*,*) 'vlspltgen 597: sortie' 
541      RETURN
542      END
Note: See TracBrowser for help on using the repository browser.