source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/config/ppsrc/dyn/vlspltgen_p.f @ 224

Last change on this file since 224 was 224, checked in by ymipsl, 10 years ago
File size: 15.8 KB
Line 
1!
2! $Header$
3!
4       SUBROUTINE vlspltgen_p( q,iadv,pente_max,masse,w,pbaru,pbarv,pdt,
5     ,                                  p,pk,teta                 )
6c
7c     Auteurs:   P.Le Van, F.Hourdin, F.Forget, F.Codron 
8c
9c    ********************************************************************
10c          Shema  d'advection " pseudo amont " .
11c      + test sur humidite specifique: Q advecte< Qsat aval
12c                   (F. Codron, 10/99)
13c    ********************************************************************
14c     q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
15c
16c     pente_max facteur de limitation des pentes: 2 en general
17c                                                0 pour un schema amont
18c     pbaru,pbarv,w flux de masse en u ,v ,w
19c     pdt pas de temps
20c
21c     teta temperature potentielle, p pression aux interfaces,
22c     pk exner au milieu des couches necessaire pour calculer Qsat
23c   --------------------------------------------------------------------
24      USE parallel_lmdz
25      USE mod_hallo
26      USE Write_Field_p
27      USE VAMPIR
28      USE infotrac, ONLY : nqtot
29      IMPLICIT NONE
30
31c
32!-----------------------------------------------------------------------
33!   INCLUDE 'dimensions.h'
34!
35!   dimensions.h contient les dimensions du modele
36!   ndm est tel que iim=2**ndm
37!-----------------------------------------------------------------------
38
39      INTEGER iim,jjm,llm,ndm
40
41      PARAMETER (iim= 128,jjm=96,llm=64,ndm=1)
42
43!-----------------------------------------------------------------------
44!
45! $Header$
46!
47!
48!  ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
49!                 veillez  n'utiliser que des ! pour les commentaires
50!                 et  bien positionner les & des lignes de continuation
51!                 (les placer en colonne 6 et en colonne 73)
52!
53!
54!-----------------------------------------------------------------------
55!   INCLUDE 'paramet.h'
56
57      INTEGER  iip1,iip2,iip3,jjp1,llmp1,llmp2,llmm1
58      INTEGER  kftd,ip1jm,ip1jmp1,ip1jmi1,ijp1llm
59      INTEGER  ijmllm,mvar
60      INTEGER jcfil,jcfllm
61
62      PARAMETER( iip1= iim+1,iip2=iim+2,iip3=iim+3                       &
63     &    ,jjp1=jjm+1-1/jjm)
64      PARAMETER( llmp1 = llm+1,  llmp2 = llm+2, llmm1 = llm-1 )
65      PARAMETER( kftd  = iim/2 -ndm )
66      PARAMETER( ip1jm  = iip1*jjm,  ip1jmp1= iip1*jjp1 )
67      PARAMETER( ip1jmi1= ip1jm - iip1 )
68      PARAMETER( ijp1llm= ip1jmp1 * llm, ijmllm= ip1jm * llm )
69      PARAMETER( mvar= ip1jmp1*( 2*llm+1) + ijmllm )
70      PARAMETER( jcfil=jjm/2+5, jcfllm=jcfil*llm )
71
72!-----------------------------------------------------------------------
73!
74! $Id: logic.h 1520 2011-05-23 11:37:09Z emillour $
75!
76!
77! NB: keep items of different kinds in seperate common blocs to avoid
78!     "misaligned commons" issues
79!-----------------------------------------------------------------------
80! INCLUDE 'logic.h'
81
82      COMMON/logicl/ purmats,forward,leapf,apphys,                      &
83     &  statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus       &
84     &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
85     &  ,ok_limit,ok_etat0,hybrid                                       &
86     &  ,moyzon_mu,moyzon_ch
87
88      COMMON/logici/ iflag_phys,iflag_trac
89     
90      LOGICAL purmats,forward,leapf,apphys,statcl,conser,               &
91     & apdiss,apdelq,saison,ecripar,fxyhypb,ysinus                      &
92     &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
93     &  ,ok_limit,ok_etat0
94      logical hybrid ! vertical coordinate is hybrid if true (sigma otherwise)
95                     ! (only used if disvert_type==2)
96      logical moyzon_mu,moyzon_ch ! used for zonal averages in Titan
97
98      integer iflag_phys,iflag_trac
99!$OMP THREADPRIVATE(/logicl/)
100!$OMP THREADPRIVATE(/logici/)
101!-----------------------------------------------------------------------
102!
103! $Id: comvert.h 1654 2012-09-24 15:07:18Z aslmd $
104!
105!-----------------------------------------------------------------------
106!   INCLUDE 'comvert.h'
107
108      COMMON/comvertr/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm),     &
109     &               pa,preff,nivsigs(llm),nivsig(llm+1),               &
110     &               aps(llm),bps(llm),scaleheight,pseudoalt(llm)
111
112      common/comverti/disvert_type, pressure_exner
113
114      real ap     ! hybrid pressure contribution at interlayers
115      real bp     ! hybrid sigma contribution at interlayer
116      real presnivs ! (reference) pressure at mid-layers
117      real dpres
118      real pa     ! reference pressure (Pa) at which hybrid coordinates
119                  ! become purely pressure
120      real preff  ! reference surface pressure (Pa)
121      real nivsigs
122      real nivsig
123      real aps    ! hybrid pressure contribution at mid-layers
124      real bps    ! hybrid sigma contribution at mid-layers
125      real scaleheight ! atmospheric (reference) scale height (km)
126      real pseudoalt ! pseudo-altitude of model levels (km), based on presnivs(),
127                     ! preff and scaleheight
128
129      integer disvert_type ! type of vertical discretization:
130                           ! 1: Earth (default for planet_type==earth),
131                           !     automatic generation
132                           ! 2: Planets (default for planet_type!=earth),
133                           !     using 'z2sig.def' (or 'esasig.def) file
134
135      logical pressure_exner
136!     compute pressure inside layers using Exner function, else use mean
137!     of pressure values at interfaces
138
139 !-----------------------------------------------------------------------
140!
141! $Id: comconst.h 1437 2010-09-30 08:29:10Z emillour $
142!
143!-----------------------------------------------------------------------
144! INCLUDE comconst.h
145
146      COMMON/comconsti/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl,          &
147     &                 iflag_top_bound,mode_top_bound
148      COMMON/comconstr/dtvr,daysec,                                     &
149     & pi,dtphys,dtdiss,rad,r,kappa,cotot,unsim,g,omeg                  &
150     & ,dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta        &
151     & ,dissip_pupstart  ,tau_top_bound,                                &
152     & daylen,molmass, ihf
153      COMMON/cpdetvenus/cpp,nu_venus,t0_venus
154
155      INTEGER im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl
156      REAL dtvr ! dynamical time step (in s)
157      REAL daysec !length (in s) of a standard day
158      REAL pi    ! something like 3.14159....
159      REAL dtphys ! (s) time step for the physics
160      REAL dtdiss ! (s) time step for the dissipation
161      REAL rad ! (m) radius of the planet
162      REAL r ! Reduced Gas constant r=R/mu
163             ! with R=8.31.. J.K-1.mol-1, mu: mol mass of atmosphere (kg/mol)
164      REAL cpp   ! Cp
165      REAL kappa ! kappa=R/Cp
166      REAL cotot
167      REAL unsim ! = 1./iim
168      REAL g ! (m/s2) gravity
169      REAL omeg ! (rad/s) rotation rate of the planet
170! Dissipation factors, for Earth model:
171      REAL dissip_factz,dissip_zref !dissip_deltaz
172! Dissipation factors, for other planets:
173      REAL dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta
174      REAL dissip_pupstart
175      INTEGER iflag_top_bound,mode_top_bound
176      REAL tau_top_bound
177      REAL daylen ! length of solar day, in 'standard' day length
178      REAL molmass ! (g/mol) molar mass of the atmosphere
179
180      REAL nu_venus,t0_venus ! coeffs needed for Cp(T), Venus atmosphere
181      REAL ihf  ! (W/m2) intrinsic heat flux for giant planets
182
183
184!-----------------------------------------------------------------------
185
186c
187c   Arguments:
188c   ----------
189      INTEGER iadv(nqtot)
190      REAL masse(ip1jmp1,llm),pente_max
191      REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
192      REAL q(ip1jmp1,llm,nqtot)
193      REAL w(ip1jmp1,llm),pdt
194      REAL p(ip1jmp1,llmp1),teta(ip1jmp1,llm),pk(ip1jmp1,llm)
195c
196c      Local 
197c   ---------
198c
199      INTEGER ij,l
200c
201      REAL,SAVE :: qsat(ip1jmp1,llm)
202      REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: zm
203      REAL,SAVE :: mu(ip1jmp1,llm)
204      REAL,SAVE :: mv(ip1jm,llm)
205      REAL,SAVE :: mw(ip1jmp1,llm+1)
206      REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: zq
207      REAL zzpbar, zzw
208
209      REAL qmin,qmax
210      DATA qmin,qmax/0.,1.e33/
211
212c--pour rapport de melange saturant--
213
214      REAL rtt,retv,r2es,r3les,r3ies,r4les,r4ies,play
215      REAL ptarg,pdelarg,foeew,zdelta
216      REAL tempe(ip1jmp1)
217      INTEGER ijb,ije,iq
218      LOGICAL, SAVE :: firstcall=.TRUE.
219!$OMP THREADPRIVATE(firstcall)
220      type(request) :: MyRequest1
221      type(request) :: MyRequest2
222
223c    fonction psat(T)
224
225       FOEEW ( PTARG,PDELARG ) = EXP (
226     *          (R3LES*(1.-PDELARG)+R3IES*PDELARG) * (PTARG-RTT)
227     * / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG)) )
228
229        r2es  = 380.11733 
230        r3les = 17.269
231        r3ies = 21.875
232        r4les = 35.86
233        r4ies = 7.66
234        retv = 0.6077667
235        rtt  = 273.16
236
237c Allocate variables depending on dynamic variable nqtot
238
239         IF (firstcall) THEN
240            firstcall=.FALSE.
241!$OMP MASTER
242            ALLOCATE(zm(ip1jmp1,llm,nqtot))
243            ALLOCATE(zq(ip1jmp1,llm,nqtot))
244!$OMP END MASTER
245!$OMP BARRIER
246         END IF
247c-- Calcul de Qsat en chaque point
248c-- approximation: au milieu des couches play(l)=(p(l)+p(l+1))/2
249c   pour eviter une exponentielle.
250
251      call SetTag(MyRequest1,100)
252      call SetTag(MyRequest2,101)
253
254       
255        ijb=ij_begin-iip1
256        ije=ij_end+iip1
257        if (pole_nord) ijb=ij_begin
258        if (pole_sud) ije=ij_end
259       
260c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 
261        DO l = 1, llm
262         DO ij = ijb, ije
263          tempe(ij) = teta(ij,l) * pk(ij,l) /cpp
264         ENDDO
265         DO ij = ijb, ije
266          zdelta = MAX( 0., SIGN(1., rtt - tempe(ij)) )
267          play   = 0.5*(p(ij,l)+p(ij,l+1))
268          qsat(ij,l) = MIN(0.5, r2es* FOEEW(tempe(ij),zdelta) / play )
269          qsat(ij,l) = qsat(ij,l) / ( 1. - retv * qsat(ij,l) )
270         ENDDO
271        ENDDO
272c$OMP END DO NOWAIT
273c      PRINT*,'Debut vlsplt version debug sans vlyqs'
274
275        zzpbar = 0.5 * pdt
276        zzw    = pdt
277
278      ijb=ij_begin
279      ije=ij_end
280      if (pole_nord) ijb=ijb+iip1
281      if (pole_sud)  ije=ije-iip1
282
283c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
284      DO l=1,llm
285        DO ij = ijb,ije
286            mu(ij,l)=pbaru(ij,l) * zzpbar
287         ENDDO
288      ENDDO
289c$OMP END DO NOWAIT
290     
291      ijb=ij_begin-iip1
292      ije=ij_end
293      if (pole_nord) ijb=ij_begin
294      if (pole_sud)  ije=ij_end-iip1
295
296c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
297      DO l=1,llm
298         DO ij=ijb,ije
299            mv(ij,l)=pbarv(ij,l) * zzpbar
300         ENDDO
301      ENDDO
302c$OMP END DO NOWAIT
303
304      ijb=ij_begin
305      ije=ij_end
306
307c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
308      DO l=1,llm
309         DO ij=ijb,ije
310            mw(ij,l)=w(ij,l) * zzw
311         ENDDO
312      ENDDO
313c$OMP END DO NOWAIT
314
315c$OMP MASTER
316      DO ij=ijb,ije
317         mw(ij,llm+1)=0.
318      ENDDO
319c$OMP END MASTER
320
321c      CALL SCOPY(ijp1llm,q,1,zq,1)
322c      CALL SCOPY(ijp1llm,masse,1,zm,1)
323
324       ijb=ij_begin
325       ije=ij_end
326
327      DO iq=1,nqtot
328c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 
329        DO l=1,llm
330          zq(ijb:ije,l,iq)=q(ijb:ije,l,iq)
331          zm(ijb:ije,l,iq)=masse(ijb:ije,l)
332        ENDDO
333c$OMP END DO NOWAIT
334      ENDDO
335
336
337c$OMP BARRIER           
338      DO iq=1,nqtot
339
340        if(iadv(iq) == 0) then
341       
342          cycle
343       
344        else if (iadv(iq)==10) then
345
346          call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
347     &               ij_begin,ij_end)
348
349c$OMP MASTER
350          call VTb(VTHallo)
351c$OMP END MASTER
352          call Register_Hallo(zq(1,1,iq),ip1jmp1,llm,2,2,2,2,MyRequest1)
353          call Register_Hallo(zm(1,1,iq),ip1jmp1,llm,1,1,1,1,MyRequest1)
354
355c$OMP MASTER
356          call VTe(VTHallo)
357c$OMP END MASTER
358        else if (iadv(iq)==14) then
359
360
361          call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat,
362     &                 ij_begin,ij_end)
363
364c$OMP MASTER
365          call VTb(VTHallo)
366c$OMP END MASTER
367
368          call Register_Hallo(zq(1,1,iq),ip1jmp1,llm,2,2,2,2,MyRequest1)
369          call Register_Hallo(zm(1,1,iq),ip1jmp1,llm,1,1,1,1,MyRequest1)
370
371c$OMP MASTER
372          call VTe(VTHallo)
373c$OMP END MASTER 
374        else
375       
376          stop 'vlspltgen_p : schema non parallelise'
377     
378        endif
379     
380      enddo
381     
382     
383c$OMP BARRIER     
384c$OMP MASTER     
385      call VTb(VTHallo)
386c$OMP END MASTER
387
388      call SendRequest(MyRequest1)
389
390c$OMP MASTER
391      call VTe(VTHallo)
392c$OMP END MASTER       
393c$OMP BARRIER
394      do iq=1,nqtot
395
396        if(iadv(iq) == 0) then
397       
398          cycle
399       
400        else if (iadv(iq)==10) then
401
402        else if (iadv(iq)==14) then
403        else
404       
405          stop 'vlspltgen_p : schema non parallelise'
406     
407        endif
408     
409      enddo
410c$OMP BARRIER     
411c$OMP MASTER
412      call VTb(VTHallo)
413c$OMP END MASTER
414
415!      call WaitRecvRequest(MyRequest1)
416!      call WaitSendRequest(MyRequest1)
417c$OMP BARRIER
418       call WaitRequest(MyRequest1)
419
420
421c$OMP MASTER
422      call VTe(VTHallo)
423c$OMP END MASTER
424c$OMP BARRIER
425 
426      do iq=1,nqtot
427
428        if(iadv(iq) == 0) then
429       
430          cycle
431       
432        else if (iadv(iq)==10) then
433       
434          call vly_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv)
435 
436        else if (iadv(iq)==14) then
437     
438          call vlyqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv,qsat)
439 
440        else
441       
442          stop 'vlspltgen_p : schema non parallelise'
443     
444        endif
445       
446       enddo
447
448
449      do iq=1,nqtot
450
451        if(iadv(iq) == 0) then
452         
453          cycle
454       
455        else if (iadv(iq)==10 .or. iadv(iq)==14 ) then
456
457c$OMP BARRIER       
458          call vlz_p(zq(1,1,iq),pente_max,zm(1,1,iq),mw,
459     &               ij_begin,ij_end)
460c$OMP BARRIER
461
462c$OMP MASTER
463          call VTb(VTHallo)
464c$OMP END MASTER
465
466          call Register_Hallo(zq(1,1,iq),ip1jmp1,llm,2,2,2,2,MyRequest2)
467          call Register_Hallo(zm(1,1,iq),ip1jmp1,llm,1,1,1,1,MyRequest2)
468
469c$OMP MASTER
470          call VTe(VTHallo)
471c$OMP END MASTER       
472c$OMP BARRIER
473        else
474       
475          stop 'vlspltgen_p : schema non parallelise'
476     
477        endif
478     
479      enddo
480c$OMP BARRIER     
481
482c$OMP MASTER       
483      call VTb(VTHallo)
484c$OMP END MASTER
485
486      call SendRequest(MyRequest2)
487
488c$OMP MASTER
489      call VTe(VTHallo)
490c$OMP END MASTER       
491
492c$OMP BARRIER
493      do iq=1,nqtot
494
495        if(iadv(iq) == 0) then
496         
497          cycle
498       
499        else if (iadv(iq)==10 .or. iadv(iq)==14 ) then
500c$OMP BARRIER       
501
502
503c$OMP BARRIER       
504        else
505       
506          stop 'vlspltgen_p : schema non parallelise'
507     
508        endif
509     
510      enddo
511
512c$OMP BARRIER
513c$OMP MASTER
514      call VTb(VTHallo)
515c$OMP END MASTER
516
517!      call WaitRecvRequest(MyRequest2)
518!      call WaitSendRequest(MyRequest2)
519c$OMP BARRIER
520       CALL WaitRequest(MyRequest2)
521
522c$OMP MASTER
523      call VTe(VTHallo)
524c$OMP END MASTER
525c$OMP BARRIER
526
527
528      do iq=1,nqtot
529
530        if(iadv(iq) == 0) then
531       
532          cycle
533       
534        else if (iadv(iq)==10) then
535       
536          call vly_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv)
537 
538        else if (iadv(iq)==14) then
539     
540          call vlyqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv,qsat)
541 
542        else
543       
544          stop 'vlspltgen_p : schema non parallelise'
545     
546        endif
547       
548       enddo
549
550      do iq=1,nqtot
551
552        if(iadv(iq) == 0) then
553         
554          cycle
555       
556        else if (iadv(iq)==10) then
557       
558          call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
559     &               ij_begin,ij_end)
560 
561        else if (iadv(iq)==14) then
562     
563          call vlxqs_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,qsat,
564     &                 ij_begin,ij_end)
565 
566        else
567       
568          stop 'vlspltgen_p : schema non parallelise'
569     
570        endif
571       
572       enddo
573
574     
575      ijb=ij_begin
576      ije=ij_end
577c$OMP BARRIER     
578
579
580      DO iq=1,nqtot
581
582c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
583        DO l=1,llm
584           DO ij=ijb,ije
585c             print *,'zq-->',ij,l,iq,zq(ij,l,iq)
586c            print *,'q-->',ij,l,iq,q(ij,l,iq)
587             q(ij,l,iq)=zq(ij,l,iq)
588           ENDDO
589        ENDDO
590c$OMP END DO NOWAIT         
591
592c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
593        DO l=1,llm
594           DO ij=ijb,ije-iip1+1,iip1
595              q(ij+iim,l,iq)=q(ij,l,iq)
596           ENDDO
597        ENDDO
598c$OMP END DO NOWAIT 
599
600      ENDDO
601       
602     
603c$OMP BARRIER
604
605cc$OMP MASTER     
606c      call WaitSendRequest(MyRequest1) 
607c      call WaitSendRequest(MyRequest2)
608cc$OMP END MASTER
609cc$OMP BARRIER
610
611      RETURN
612      END
Note: See TracBrowser for help on using the repository browser.