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

Last change on this file since 224 was 224, checked in by ymipsl, 10 years ago
File size: 15.7 KB
Line 
1!
2! $Header$
3!
4      SUBROUTINE advect_new_p(ucov,vcov,teta,w,massebx,masseby,
5     &                        du,dv,dteta)
6      USE parallel_lmdz
7      USE write_field_p
8      IMPLICIT NONE
9c=======================================================================
10c
11c   Auteurs:  P. Le Van , Fr. Hourdin  .
12c   -------
13c
14c   Objet:
15c   ------
16c
17c   *************************************************************
18c   .... calcul des termes d'advection vertic.pour u,v,teta,q ...
19c   *************************************************************
20c        ces termes sont ajoutes a du,dv,dteta et dq .
21c  Modif F.Forget 03/94 : on retire q de advect
22c
23c=======================================================================
24c-----------------------------------------------------------------------
25c   Declarations:
26c   -------------
27
28!-----------------------------------------------------------------------
29!   INCLUDE 'dimensions.h'
30!
31!   dimensions.h contient les dimensions du modele
32!   ndm est tel que iim=2**ndm
33!-----------------------------------------------------------------------
34
35      INTEGER iim,jjm,llm,ndm
36
37      PARAMETER (iim= 128,jjm=96,llm=64,ndm=1)
38
39!-----------------------------------------------------------------------
40!
41! $Header$
42!
43!
44!  ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
45!                 veillez  n'utiliser que des ! pour les commentaires
46!                 et  bien positionner les & des lignes de continuation
47!                 (les placer en colonne 6 et en colonne 73)
48!
49!
50!-----------------------------------------------------------------------
51!   INCLUDE 'paramet.h'
52
53      INTEGER  iip1,iip2,iip3,jjp1,llmp1,llmp2,llmm1
54      INTEGER  kftd,ip1jm,ip1jmp1,ip1jmi1,ijp1llm
55      INTEGER  ijmllm,mvar
56      INTEGER jcfil,jcfllm
57
58      PARAMETER( iip1= iim+1,iip2=iim+2,iip3=iim+3                       &
59     &    ,jjp1=jjm+1-1/jjm)
60      PARAMETER( llmp1 = llm+1,  llmp2 = llm+2, llmm1 = llm-1 )
61      PARAMETER( kftd  = iim/2 -ndm )
62      PARAMETER( ip1jm  = iip1*jjm,  ip1jmp1= iip1*jjp1 )
63      PARAMETER( ip1jmi1= ip1jm - iip1 )
64      PARAMETER( ijp1llm= ip1jmp1 * llm, ijmllm= ip1jm * llm )
65      PARAMETER( mvar= ip1jmp1*( 2*llm+1) + ijmllm )
66      PARAMETER( jcfil=jjm/2+5, jcfllm=jcfil*llm )
67
68!-----------------------------------------------------------------------
69!
70! $Id: comconst.h 1437 2010-09-30 08:29:10Z emillour $
71!
72!-----------------------------------------------------------------------
73! INCLUDE comconst.h
74
75      COMMON/comconsti/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl,          &
76     &                 iflag_top_bound,mode_top_bound
77      COMMON/comconstr/dtvr,daysec,                                     &
78     & pi,dtphys,dtdiss,rad,r,kappa,cotot,unsim,g,omeg                  &
79     & ,dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta        &
80     & ,dissip_pupstart  ,tau_top_bound,                                &
81     & daylen,molmass, ihf
82      COMMON/cpdetvenus/cpp,nu_venus,t0_venus
83
84      INTEGER im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl
85      REAL dtvr ! dynamical time step (in s)
86      REAL daysec !length (in s) of a standard day
87      REAL pi    ! something like 3.14159....
88      REAL dtphys ! (s) time step for the physics
89      REAL dtdiss ! (s) time step for the dissipation
90      REAL rad ! (m) radius of the planet
91      REAL r ! Reduced Gas constant r=R/mu
92             ! with R=8.31.. J.K-1.mol-1, mu: mol mass of atmosphere (kg/mol)
93      REAL cpp   ! Cp
94      REAL kappa ! kappa=R/Cp
95      REAL cotot
96      REAL unsim ! = 1./iim
97      REAL g ! (m/s2) gravity
98      REAL omeg ! (rad/s) rotation rate of the planet
99! Dissipation factors, for Earth model:
100      REAL dissip_factz,dissip_zref !dissip_deltaz
101! Dissipation factors, for other planets:
102      REAL dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta
103      REAL dissip_pupstart
104      INTEGER iflag_top_bound,mode_top_bound
105      REAL tau_top_bound
106      REAL daylen ! length of solar day, in 'standard' day length
107      REAL molmass ! (g/mol) molar mass of the atmosphere
108
109      REAL nu_venus,t0_venus ! coeffs needed for Cp(T), Venus atmosphere
110      REAL ihf  ! (W/m2) intrinsic heat flux for giant planets
111
112
113!-----------------------------------------------------------------------
114!
115! $Id: comvert.h 1654 2012-09-24 15:07:18Z aslmd $
116!
117!-----------------------------------------------------------------------
118!   INCLUDE 'comvert.h'
119
120      COMMON/comvertr/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm),     &
121     &               pa,preff,nivsigs(llm),nivsig(llm+1),               &
122     &               aps(llm),bps(llm),scaleheight,pseudoalt(llm)
123
124      common/comverti/disvert_type, pressure_exner
125
126      real ap     ! hybrid pressure contribution at interlayers
127      real bp     ! hybrid sigma contribution at interlayer
128      real presnivs ! (reference) pressure at mid-layers
129      real dpres
130      real pa     ! reference pressure (Pa) at which hybrid coordinates
131                  ! become purely pressure
132      real preff  ! reference surface pressure (Pa)
133      real nivsigs
134      real nivsig
135      real aps    ! hybrid pressure contribution at mid-layers
136      real bps    ! hybrid sigma contribution at mid-layers
137      real scaleheight ! atmospheric (reference) scale height (km)
138      real pseudoalt ! pseudo-altitude of model levels (km), based on presnivs(),
139                     ! preff and scaleheight
140
141      integer disvert_type ! type of vertical discretization:
142                           ! 1: Earth (default for planet_type==earth),
143                           !     automatic generation
144                           ! 2: Planets (default for planet_type!=earth),
145                           !     using 'z2sig.def' (or 'esasig.def) file
146
147      logical pressure_exner
148!     compute pressure inside layers using Exner function, else use mean
149!     of pressure values at interfaces
150
151 !-----------------------------------------------------------------------
152!
153! $Header$
154!
155!CDK comgeom
156      COMMON/comgeom/                                                   &
157     & cu(ip1jmp1),cv(ip1jm),unscu2(ip1jmp1),unscv2(ip1jm),             &
158     & aire(ip1jmp1),airesurg(ip1jmp1),aireu(ip1jmp1),                  &
159     & airev(ip1jm),unsaire(ip1jmp1),apoln,apols,                       &
160     & unsairez(ip1jm),airuscv2(ip1jm),airvscu2(ip1jm),                 &
161     & aireij1(ip1jmp1),aireij2(ip1jmp1),aireij3(ip1jmp1),              &
162     & aireij4(ip1jmp1),alpha1(ip1jmp1),alpha2(ip1jmp1),                &
163     & alpha3(ip1jmp1),alpha4(ip1jmp1),alpha1p2(ip1jmp1),               &
164     & alpha1p4(ip1jmp1),alpha2p3(ip1jmp1),alpha3p4(ip1jmp1),           &
165     & fext(ip1jm),constang(ip1jmp1),rlatu(jjp1),rlatv(jjm),            &
166     & rlonu(iip1),rlonv(iip1),cuvsurcv(ip1jm),cvsurcuv(ip1jm),         &
167     & cvusurcu(ip1jmp1),cusurcvu(ip1jmp1),cuvscvgam1(ip1jm),           &
168     & cuvscvgam2(ip1jm),cvuscugam1(ip1jmp1),                           &
169     & cvuscugam2(ip1jmp1),cvscuvgam(ip1jm),cuscvugam(ip1jmp1),         &
170     & unsapolnga1,unsapolnga2,unsapolsga1,unsapolsga2,                 &
171     & unsair_gam1(ip1jmp1),unsair_gam2(ip1jmp1),unsairz_gam(ip1jm),    &
172     & aivscu2gam(ip1jm),aiuscv2gam(ip1jm),xprimu(iip1),xprimv(iip1)
173
174!
175        REAL                                                            &
176     & cu,cv,unscu2,unscv2,aire,airesurg,aireu,airev,unsaire,apoln     ,&
177     & apols,unsairez,airuscv2,airvscu2,aireij1,aireij2,aireij3,aireij4,&
178     & alpha1,alpha2,alpha3,alpha4,alpha1p2,alpha1p4,alpha2p3,alpha3p4 ,&
179     & fext,constang,rlatu,rlatv,rlonu,rlonv,cuvscvgam1,cuvscvgam2     ,&
180     & cvuscugam1,cvuscugam2,cvscuvgam,cuscvugam,unsapolnga1,unsapolnga2&
181     & ,unsapolsga1,unsapolsga2,unsair_gam1,unsair_gam2,unsairz_gam    ,&
182     & aivscu2gam ,aiuscv2gam,cuvsurcv,cvsurcuv,cvusurcu,cusurcvu,xprimu&
183     & , xprimv
184!
185!
186! $Id: logic.h 1520 2011-05-23 11:37:09Z emillour $
187!
188!
189! NB: keep items of different kinds in seperate common blocs to avoid
190!     "misaligned commons" issues
191!-----------------------------------------------------------------------
192! INCLUDE 'logic.h'
193
194      COMMON/logicl/ purmats,forward,leapf,apphys,                      &
195     &  statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus       &
196     &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
197     &  ,ok_limit,ok_etat0,hybrid                                       &
198     &  ,moyzon_mu,moyzon_ch
199
200      COMMON/logici/ iflag_phys,iflag_trac
201     
202      LOGICAL purmats,forward,leapf,apphys,statcl,conser,               &
203     & apdiss,apdelq,saison,ecripar,fxyhypb,ysinus                      &
204     &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
205     &  ,ok_limit,ok_etat0
206      logical hybrid ! vertical coordinate is hybrid if true (sigma otherwise)
207                     ! (only used if disvert_type==2)
208      logical moyzon_mu,moyzon_ch ! used for zonal averages in Titan
209
210      integer iflag_phys,iflag_trac
211!$OMP THREADPRIVATE(/logicl/)
212!$OMP THREADPRIVATE(/logici/)
213!-----------------------------------------------------------------------
214!
215! $Id: ener.h 1447 2010-10-22 16:18:27Z jghattas $
216!
217!  ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
218!                 veillez à n'utiliser que des ! pour les commentaires
219!                 et à bien positionner les & des lignes de continuation
220!                 (les placer en colonne 6 et en colonne 73)
221!
222! INCLUDE 'ener.h'
223
224      COMMON/ener/ang0,etot0,ptot0,ztot0,stot0,                         &
225     &            ang,etot,ptot,ztot,stot,rmsdpdt ,                     &
226     &            rmsv,gtot(llmm1)
227
228      REAL ang0,etot0,ptot0,ztot0,stot0,                                &
229     &     ang,etot,ptot,ztot,stot,rmsdpdt,rmsv,gtot
230
231!-----------------------------------------------------------------------
232
233c   Arguments:
234c   ----------
235
236      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm),teta(ip1jmp1,llm)
237      REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm),w(ip1jmp1,llm)
238      REAL dv(ip1jm,llm),du(ip1jmp1,llm),dteta(ip1jmp1,llm)
239      REAL,SAVE :: dv1(ip1jm,llm),du1(ip1jmp1,llm),dteta1(ip1jmp1,llm)
240      REAL,SAVE :: dv2(ip1jm,llm),du2(ip1jmp1,llm),dteta2(ip1jmp1,llm)
241c   Local:
242c   ------
243
244      REAL,SAVE :: uav(ip1jmp1,llm),vav(ip1jm,llm)
245      REAL wsur2(ip1jmp1)
246      REAL unsaire2(ip1jmp1), ge(ip1jmp1)
247      REAL deuxjour, ww, gt, uu, vv
248
249      INTEGER  ij,l,ijb,ije
250
251      EXTERNAL  SSUM
252      REAL      SSUM
253
254c-----------------------------------------------------------------------
255c   2. Calculs preliminaires:
256c   -------------------------
257
258      IF (conser)  THEN
259         deuxjour = 2. * daysec
260
261         DO   1  ij   = 1, ip1jmp1
262         unsaire2(ij) = unsaire(ij) * unsaire(ij)
263   1     CONTINUE
264      END IF
265
266
267c------------------  -yy ----------------------------------------------
268c   .  Calcul de     u
269
270c$OMP MASTER
271      ijb=ij_begin
272      ije=ij_end
273      if (pole_nord) ijb=ijb+iip1
274      if (pole_sud)  ije=ije-iip1
275
276      DO ij=ijb,ije
277        du2(ij,1)=0.
278      ENDDO
279     
280      ijb=ij_begin
281      ije=ij_end
282      if (pole_sud)  ije=ij_end-iip1
283     
284      DO ij=ijb,ije
285        dv2(ij,1)=0.
286      ENDDO
287     
288      ijb=ij_begin
289      ije=ij_end
290
291      DO ij=ijb,ije
292        dteta2(ij,1)=0.
293      ENDDO
294c$OMP END MASTER
295
296 
297c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
298      DO  l=1,llm
299         
300         ijb=ij_begin
301         ije=ij_end
302         if (pole_nord) ijb=ijb+iip1
303         if (pole_sud)  ije=ije-iip1
304         
305c         DO    ij     = iip2, ip1jmp1
306c            uav(ij,l) = 0.25 * ( ucov(ij,l) + ucov(ij-iip1,l) )
307c         ENDDO
308
309c         DO    ij     = iip2, ip1jm
310c            uav(ij,l) = uav(ij,l) + uav(ij+iip1,l)
311c         ENDDO
312         
313         DO    ij     = ijb, ije
314                 
315           uav(ij,l)=0.25*(ucov(ij,l)+ucov(ij-iip1,l))
316     .               +0.25*(ucov(ij+iip1,l)+ucov(ij,l))
317         ENDDO
318         
319         if (pole_nord) then
320           DO      ij         = 1, iip1
321              uav(ij      ,l) = 0.
322           ENDDO
323         endif
324         
325         if (pole_sud) then
326           DO      ij         = 1, iip1
327              uav(ip1jm+ij,l) = 0.
328           ENDDO
329         endif
330         
331      ENDDO
332c$OMP END DO     
333c      call write_field3d_p('uav',reshape(uav,(/iip1,jjp1,llm/)))
334     
335c------------------  -xx ----------------------------------------------
336c   .  Calcul de     v
337     
338      ijb=ij_begin
339      ije=ij_end
340      if (pole_sud)  ije=ij_end-iip1
341
342c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
343      DO  l=1,llm
344         
345         DO    ij   = ijb+1, ije
346           vav(ij,l) = 0.25 * ( vcov(ij,l) + vcov(ij-1,l) )
347         ENDDO
348         
349         DO    ij   = ijb,ije,iip1
350          vav(ij,l) = vav(ij+iim,l)
351         ENDDO
352         
353         
354         DO    ij   = ijb, ije-1
355          vav(ij,l) = vav(ij,l) + vav(ij+1,l)
356         ENDDO
357         
358         DO    ij       = ijb, ije, iip1
359          vav(ij+iim,l) = vav(ij,l)
360         ENDDO
361         
362      ENDDO
363c$OMP END DO
364c       call write_field3d_p('vav',reshape(vav,(/iip1,jjm,llm/)))
365
366c-----------------------------------------------------------------------
367c$OMP BARRIER
368
369c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
370      DO 20 l = 1, llmm1
371
372
373c       ......   calcul de  - w/2.    au niveau  l+1   .......
374      ijb=ij_begin
375      ije=ij_end+iip1
376      if (pole_sud)  ije=ij_end
377     
378      DO 5   ij   = ijb, ije
379      wsur2( ij ) = - 0.5 * w( ij,l+1 )
380   5  CONTINUE
381
382
383c     .....................     calcul pour  du     ..................
384     
385      ijb=ij_begin
386      ije=ij_end
387      if (pole_nord) ijb=ijb+iip1
388      if (pole_sud)  ije=ije-iip1
389         
390      DO 6 ij = ijb ,ije-1
391      ww        = wsur2 (  ij  )     + wsur2( ij+1 ) 
392      uu        = 0.5 * ( ucov(ij,l) + ucov(ij,l+1) )
393      du1(ij,l)  =  ww * ( uu - uav(ij, l ) )/massebx(ij, l )
394      du2(ij,l+1)=  ww * ( uu - uav(ij,l+1) )/massebx(ij,l+1)
395   6  CONTINUE
396
397c     .................    calcul pour   dv      .....................
398      ijb=ij_begin
399      ije=ij_end
400      if (pole_sud)  ije=ij_end-iip1
401     
402      DO 8 ij = ijb, ije
403      ww        = wsur2( ij+iip1 )   + wsur2( ij )
404      vv        = 0.5 * ( vcov(ij,l) + vcov(ij,l+1) )
405      dv1(ij,l)  =  ww * (vv - vav(ij, l ) )/masseby(ij, l )
406      dv2(ij,l+1)=  ww * (vv - vav(ij,l+1) )/masseby(ij,l+1)
407   8  CONTINUE
408
409c
410
411c     ............................................................
412c     ...............    calcul pour   dh      ...................
413c     ............................................................
414
415c                       ---z
416c       calcul de  - d( teta  * w )      qu'on ajoute a   dh
417c                   ...............
418        ijb=ij_begin
419        ije=ij_end
420       
421        DO 15 ij = ijb, ije
422         ww            = wsur2(ij) * (teta(ij,l) + teta(ij,l+1) )
423         dteta1(ij, l ) =   ww
424         dteta2(ij,l+1) =   ww
425  15    CONTINUE
426
427c ym ---> conser a voir plus tard
428
429c      IF( conser)  THEN
430c       
431c        DO 17 ij = 1,ip1jmp1
432c        ge(ij)   = wsur2(ij) * wsur2(ij) * unsaire2(ij)
433c  17    CONTINUE
434c        gt       = SSUM( ip1jmp1,ge,1 )
435c        gtot(l)  = deuxjour * SQRT( gt/ip1jmp1 )
436c      END IF
437
438  20  CONTINUE
439c$OMP END DO
440
441      ijb=ij_begin
442      ije=ij_end
443      if (pole_nord) ijb=ijb+iip1
444      if (pole_sud)  ije=ije-iip1
445     
446c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
447      DO l=1,llm
448        DO ij=ijb,ije-1
449          du(ij,l)=du(ij,l)+du2(ij,l)-du1(ij,l)
450        ENDDO
451
452        DO   ij   = ijb+iip1-1, ije, iip1
453         du( ij, l  ) = du( ij -iim, l  )
454        ENDDO 
455      ENDDO
456c$OMP END DO NOWAIT
457     
458      ijb=ij_begin
459      ije=ij_end
460      if (pole_sud)  ije=ij_end-iip1
461
462c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
463      DO l=1,llm
464        DO ij=ijb,ije
465          dv(ij,l)=dv(ij,l)+dv2(ij,l)-dv1(ij,l)
466        ENDDO
467      ENDDO
468c$OMP END DO NOWAIT     
469      ijb=ij_begin
470      ije=ij_end
471
472c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
473      DO l=1,llm
474        DO ij=ijb,ije
475          dteta(ij,l)=dteta(ij,l)+dteta2(ij,l)-dteta1(ij,l)
476        ENDDO
477      ENDDO
478c$OMP END DO NOWAIT     
479
480      RETURN
481      END
Note: See TracBrowser for help on using the repository browser.