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

Last change on this file since 5087 was 4996, checked in by evignon, 12 months ago

ajout d'un flag pour le calcul de qsat dans la condtion de "francis"
pour l'advection de l'humidite (q<qsat_aval). En activant ce flag,
on calcule qsat /liquide quelque soit la temperature et on peut donc
ainsi autoriser l'advection de sursaturations / glace à T<0oC.
  • 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.6 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          Schema  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      USE logic_mod, ONLY: adv_qsat_liq
35      IMPLICIT NONE
36
37c
38      include "dimensions.h"
39      include "paramet.h"
40
41c
42c   Arguments:
43c   ----------
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)
66      INTEGER ijb,ije,iq,iq2,ifils
67      LOGICAL, SAVE :: firstcall=.TRUE.
68!$OMP THREADPRIVATE(firstcall)
69      type(request),SAVE :: MyRequest1
70!$OMP THREADPRIVATE(MyRequest1)
71      type(request),SAVE :: MyRequest2
72!$OMP THREADPRIVATE(MyRequest2)
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       
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       
105c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
106        DO l = 1, llm
107         DO ij = ijb, ije
108          tempe(ij) = teta(ij,l) * pk(ij,l) /cpp
109         ENDDO
110         DO ij = ijb, ije
111          IF (adv_qsat_liq) THEN
112             zdelta = 0.
113          ELSE
114             zdelta = MAX( 0., SIGN(1., rtt - tempe(ij)) )
115          ENDIF
116          play   = 0.5*(p(ij,l)+p(ij,l+1))
117          qsat(ij,l) = MIN(0.5, r2es* FOEEW(tempe(ij),zdelta) / play )
118          qsat(ij,l) = qsat(ij,l) / ( 1. - retv * qsat(ij,l) )
119         ENDDO
120        ENDDO
121c$OMP END DO NOWAIT
122c      PRINT*,'Debut vlsplt version debug sans vlyqs'
123
124        zzpbar = 0.5 * pdt
125        zzw    = pdt
126
127      ijb=ij_begin
128      ije=ij_end
129      if (pole_nord) ijb=ijb+iip1
130      if (pole_sud)  ije=ije-iip1
131
132c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
133      DO l=1,llm
134        DO ij = ijb,ije
135            mu(ij,l)=pbaru(ij,l) * zzpbar
136         ENDDO
137      ENDDO
138c$OMP END DO NOWAIT
139     
140      ijb=ij_begin-iip1
141      ije=ij_end
142      if (pole_nord) ijb=ij_begin
143      if (pole_sud)  ije=ij_end-iip1
144
145c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
146      DO l=1,llm
147         DO ij=ijb,ije
148            mv(ij,l)=pbarv(ij,l) * zzpbar
149         ENDDO
150      ENDDO
151c$OMP END DO NOWAIT
152
153      ijb=ij_begin
154      ije=ij_end
155
156      DO iq=1,nqtot
157c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
158      DO l=1,llm
159         DO ij=ijb,ije
160            mw(ij,l,iq)=w(ij,l) * zzw
161         ENDDO
162      ENDDO
163c$OMP END DO NOWAIT
164      ENDDO
165
166      DO iq=1,nqtot 
167c$OMP MASTER
168      DO ij=ijb,ije
169         mw(ij,llm+1,iq)=0.
170      ENDDO
171c$OMP END MASTER
172      ENDDO
173
174c      CALL SCOPY(ijp1llm,q,1,zq,1)
175c      CALL SCOPY(ijp1llm,masse,1,zm,1)
176
177       ijb=ij_begin
178       ije=ij_end
179
180      DO iq=1,nqtot       
181c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
182        DO l=1,llm
183          zq(ijb:ije,l,iq)=q(ijb:ije,l,iq)
184          zm(ijb:ije,l,iq)=masse(ijb:ije,l)
185        ENDDO
186c$OMP END DO NOWAIT
187      ENDDO
188
189#ifdef DEBUG_IO     
190       CALL WriteField_u('mu',mu)
191       CALL WriteField_v('mv',mv)
192       CALL WriteField_u('mw',mw)
193       CALL WriteField_u('qsat',qsat)
194#endif
195
196      ! verif temporaire
197      ijb=ij_begin
198      ije=ij_end 
199      call check_isotopes(zq,ijb,ije,'vlspltgen_loc 191')
200
201c$OMP BARRIER           
202      DO iq=1,nqtot
203        ! CRisi: on ne boucle que sur les pères= ceux qui sont transportés directement par l'air
204        IF(tracers(iq)%parent /= 'air') CYCLE
205        !write(*,*) 'vlspltgen 192: iq,iadv=',iq,tracers(iq)%iadv
206#ifdef DEBUG_IO   
207        CALL WriteField_u('zq',zq(:,:,iq))
208        CALL WriteField_u('zm',zm(:,:,iq))
209#endif
210        SELECT CASE(tracers(iq)%iadv)
211          CASE(0); CYCLE
212          CASE(10)
213#ifdef _ADV_HALO       
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,
217     &                     ij_begin,ij_begin+2*iip1-1,iq)
218          call vlx_loc(zq,pente_max,zm,mu,
219     &               ij_end-2*iip1+1,ij_end,iq)
220#else
221          call vlx_loc(zq,pente_max,zm,mu,
222     &                     ij_begin,ij_end,iq)
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)
230! CRisi
231          do ifils=1,tracers(iq)%nqDescen
232            iq2=tracers(iq)%iqDescen(ifils)
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
236
237c$OMP MASTER
238          call VTe(VTHallo)
239c$OMP END MASTER
240          CASE(14)
241#ifdef _ADV_HALO           
242          call vlxqs_loc(zq,pente_max,zm,mu,
243     &                   qsat,ij_begin,ij_begin+2*iip1-1,iq)
244          call vlxqs_loc(zq,pente_max,zm,mu,
245     &                   qsat,ij_end-2*iip1+1,ij_end,iq)
246#else
247          call vlxqs_loc(zq,pente_max,zm,mu,
248     &                   qsat,ij_begin,ij_end,iq)
249#endif
250
251c$OMP MASTER
252          call VTb(VTHallo)
253c$OMP END MASTER
254
255          call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest1)
256          call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest1)
257          do ifils=1,tracers(iq)%nqDescen
258            iq2=tracers(iq)%iqDescen(ifils)
259            call Register_Hallo_u(zq(:,:,iq2),llm,2,2,2,2,MyRequest1)
260            call Register_Hallo_u(zm(:,:,iq2),llm,1,1,1,1,MyRequest1)
261          enddo
262
263c$OMP MASTER
264          call VTe(VTHallo)
265c$OMP END MASTER
266          CASE DEFAULT
267             CALL abort_gcm("vlspltgen_loc","schema non parallelise",1)
268        END SELECT
269     
270      enddo !DO iq=1,nqtot
271     
272     
273c$OMP BARRIER     
274c$OMP MASTER     
275      call VTb(VTHallo)
276c$OMP END MASTER
277
278      call SendRequest(MyRequest1)
279
280c$OMP MASTER
281      call VTe(VTHallo)
282c$OMP END MASTER       
283c$OMP BARRIER
284
285      ! verif temporaire
286      ijb=ij_begin-2*iip1
287      ije=ij_end+2*iip1 
288      if (pole_nord) ijb=ij_begin
289      if (pole_sud)  ije=ij_end 
290      call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 280')
291
292      do iq=1,nqtot
293        IF(tracers(iq)%parent /= 'air') CYCLE
294        !write(*,*) 'vlspltgen 279: iq=',iq
295
296        SELECT CASE(tracers(iq)%iadv)
297          CASE(0); CYCLE
298          CASE(10)
299#ifdef _ADV_HALLO
300          call vlx_loc(zq,pente_max,zm,mu,
301     &                 ij_begin+2*iip1,ij_end-2*iip1,iq)
302#endif       
303          CASE(14)
304#ifdef _ADV_HALLO
305          call vlxqs_loc(zq,pente_max,zm,mu,
306     &                    qsat,ij_begin+2*iip1,ij_end-2*iip1,iq)
307#endif   
308          CASE DEFAULT
309          CALL abort_gcm("vlspltgen_p","schema non parallelise",1)
310        END SELECT
311     
312      enddo
313c$OMP BARRIER     
314c$OMP MASTER
315      call VTb(VTHallo)
316c$OMP END MASTER
317
318!      call WaitRecvRequest(MyRequest1)
319!      call WaitSendRequest(MyRequest1)
320c$OMP BARRIER
321       call WaitRequest(MyRequest1)
322
323
324c$OMP MASTER
325      call VTe(VTHallo)
326c$OMP END MASTER
327c$OMP BARRIER
328
329     
330      IF(isoCheck) THEN
331           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 326')
332           ijb=ij_begin-2*iip1
333           ije=ij_end+2*iip1
334           if (pole_nord) ijb=ij_begin
335           if (pole_sud)  ije=ij_end
336           call check_isotopes(zq,ijb,ije,'vlspltgen_loc 336')
337      END IF
338
339      do iq = 1, nqtot
340       IF(tracers(iq)%parent /= 'air') CYCLE
341       !write(*,*) 'vlspltgen 321: iq=',iq
342#ifdef DEBUG_IO   
343       CALL WriteField_u('zq',zq(:,:,iq))
344       CALL WriteField_u('zm',zm(:,:,iq))
345#endif
346
347        SELECT CASE(tracers(iq)%iadv)
348          CASE(0); CYCLE
349          CASE(10); call   vly_loc(zq,pente_max,zm,mv,     iq)
350          CASE(14); call vlyqs_loc(zq,pente_max,zm,mv,qsat,iq)
351          CASE DEFAULT
352          CALL abort_gcm("vlspltgen_p","schema non parallelise",1)
353        END SELECT
354       
355       enddo
356
357      call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 357')
358
359      do iq = 1, nqtot
360       IF(tracers(iq)%parent /= 'air') CYCLE
361      !write(*,*) 'vlspltgen 349: iq=',iq
362#ifdef DEBUG_IO   
363       CALL WriteField_u('zq',zq(:,:,iq))
364       CALL WriteField_u('zm',zm(:,:,iq))
365#endif
366        SELECT CASE(tracers(iq)%iadv)
367          CASE(0); CYCLE
368          CASE(10,14)
369c$OMP BARRIER       
370#ifdef _ADV_HALLO
371          call vlz_loc(zq,pente_max,zm,mw,
372     &               ij_begin,ij_begin+2*iip1-1,iq)
373          call vlz_loc(zq,pente_max,zm,mw,
374     &               ij_end-2*iip1+1,ij_end,iq)
375#else
376          call vlz_loc(zq,pente_max,zm,mw,
377     &               ij_begin,ij_end,iq)
378#endif
379c$OMP BARRIER
380
381c$OMP MASTER
382          call VTb(VTHallo)
383c$OMP END MASTER
384
385          call Register_Hallo_u(zq(:,:,iq),llm,2,2,2,2,MyRequest2)
386          call Register_Hallo_u(zm(:,:,iq),llm,1,1,1,1,MyRequest2)
387          ! CRisi
388          do ifils=1,tracers(iq)%nqDescen
389            iq2=tracers(iq)%iqDescen(ifils)
390            call Register_Hallo_u(zq(:,:,iq2),llm,2,2,2,2,MyRequest2)
391            call Register_Hallo_u(zm(:,:,iq2),llm,1,1,1,1,MyRequest2)
392          enddo     
393c$OMP MASTER
394          call VTe(VTHallo)
395c$OMP END MASTER       
396c$OMP BARRIER
397          CASE DEFAULT
398           
399            CALL abort_gcm("vlspltgen_p","schema non parallelise",1)
400         END SELECT
401     
402      enddo
403c$OMP BARRIER     
404
405c$OMP MASTER       
406      call VTb(VTHallo)
407c$OMP END MASTER
408
409      call SendRequest(MyRequest2)
410
411c$OMP MASTER
412      call VTe(VTHallo)
413c$OMP END MASTER       
414
415
416      call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 429')
417
418c$OMP BARRIER
419      do iq=1,nqtot
420        IF(tracers(iq)%parent /= 'air') CYCLE
421      !write(*,*) 'vlspltgen 409: iq=',iq
422
423        SELECT CASE(tracers(iq)%iadv)
424          CASE(0); CYCLE
425          CASE(10,14)
426c$OMP BARRIER       
427
428#ifdef _ADV_HALLO
429          call vlz_loc(zq,pente_max,zm,mw,
430     &               ij_begin+2*iip1,ij_end-2*iip1,iq)
431#endif
432
433c$OMP BARRIER       
434          CASE DEFAULT
435          CALL abort_gcm("vlspltgen_p","schema non parallelise",1)
436        END SELECT
437     
438      enddo
439      !write(*,*) 'vlspltgen_loc 476'
440
441c$OMP BARRIER
442      !write(*,*) 'vlspltgen_loc 477'
443c$OMP MASTER
444      call VTb(VTHallo)
445c$OMP END MASTER
446
447!      call WaitRecvRequest(MyRequest2)
448!      call WaitSendRequest(MyRequest2)
449c$OMP BARRIER
450       CALL WaitRequest(MyRequest2)
451
452c$OMP MASTER
453      call VTe(VTHallo)
454c$OMP END MASTER
455c$OMP BARRIER
456
457
458      !write(*,*) 'vlspltgen_loc 494'
459      call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 461')
460
461      do iq=1,nqtot
462        IF(tracers(iq)%parent /= 'air') CYCLE
463      !write(*,*) 'vlspltgen 449: iq=',iq
464#ifdef DEBUG_IO   
465       CALL WriteField_u('zq',zq(:,:,iq))
466       CALL WriteField_u('zm',zm(:,:,iq))
467#endif
468        SELECT CASE(tracers(iq)%iadv)
469          CASE(0); CYCLE
470          CASE(10); call   vly_loc(zq,pente_max,zm,mv,     iq)
471          CASE(14); call vlyqs_loc(zq,pente_max,zm,mv,qsat,iq)
472          CASE DEFAULT
473             CALL abort_gcm("vlspltgen_p","schema non parallelise",1)
474        END SELECT
475       
476       enddo !do iq=1,nqtot
477
478      call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 493')
479
480      do iq=1,nqtot
481        IF(tracers(iq)%parent /= 'air') CYCLE
482      !write(*,*) 'vlspltgen 477: iq=',iq
483#ifdef DEBUG_IO   
484       CALL WriteField_u('zq',zq(:,:,iq))
485       CALL WriteField_u('zm',zm(:,:,iq))
486#endif
487        SELECT CASE(tracers(iq)%iadv)
488          CASE(0); CYCLE
489          CASE(10); call   vlx_loc(zq,pente_max,zm,mu,
490     &               ij_begin,ij_end,iq)
491          CASE(14); call vlxqs_loc(zq,pente_max,zm,mu,
492     &                 qsat, ij_begin,ij_end,iq)
493          CASE DEFAULT
494          CALL abort_gcm("vlspltgen_p","schema non parallelise",1)
495        END SELECT
496       
497       enddo !do iq=1,nqtot
498
499      !write(*,*) 'vlspltgen 550: apres derniere serie de call vlx'
500      call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 521')
501     
502      ijb=ij_begin
503      ije=ij_end
504      !write(*,*) 'vlspltgen_loc 557'
505c$OMP BARRIER     
506
507      !write(*,*) 'vlspltgen_loc 559' 
508      DO iq=1,nqtot
509       !write(*,*) 'vlspltgen_loc 561, iq=',iq 
510#ifdef DEBUG_IO   
511       CALL WriteField_u('zq',zq(:,:,iq))
512       CALL WriteField_u('zm',zm(:,:,iq))
513#endif
514c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
515        DO l=1,llm
516           DO ij=ijb,ije
517c             print *,'zq-->',ij,l,iq,zq(ij,l,iq)
518c             print *,'q-->',ij,l,iq,q(ij,l,iq)
519             q(ij,l,iq)=zq(ij,l,iq)
520           ENDDO
521        ENDDO
522c$OMP END DO NOWAIT   
523      !write(*,*) 'vlspltgen_loc 575'     
524
525c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
526        DO l=1,llm
527           DO ij=ijb,ije-iip1+1,iip1
528              q(ij+iim,l,iq)=q(ij,l,iq)
529           ENDDO
530        ENDDO
531c$OMP END DO NOWAIT 
532      !write(*,*) 'vlspltgen_loc 583' 
533      ENDDO !DO iq=1,nqtot
534       
535      call check_isotopes(q,ij_begin,ij_end,'vlspltgen_loc 557')
536
537c$OMP BARRIER
538
539cc$OMP MASTER     
540c      call WaitSendRequest(MyRequest1)
541c      call WaitSendRequest(MyRequest2)
542cc$OMP END MASTER
543cc$OMP BARRIER
544
545      !write(*,*) 'vlspltgen 597: sortie' 
546      RETURN
547      END
Note: See TracBrowser for help on using the repository browser.