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

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