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

Last change on this file since 5183 was 4996, checked in by evignon, 6 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.