trunk/dyn3d/PVtheta.f

      SUBROUTINE PVtheta(ilon,ilev,pucov,pvcov,pteta,
     &     ztfi,zplay,zplev,
     &     nbteta,theta,PVteta)
      use dimens_m
      use paramet_m
      use comconst
      use disvert_m
      use comgeom
      use tourabs_m, only: tourabs
      IMPLICIT none

c=======================================================================
c
c   Auteur:  I. Musat
c   -------
c
c   Objet:
c   ------
c
c    *******************************************************************
c    Calcul de la vorticite potentielle PVteta sur des iso-theta selon
c    la methodologie du NCEP/NCAR :
c    1) on calcule la stabilite statique N**2=g/T*(dT/dz+g/cp) sur les
c       niveaux du modele => N2
c    2) on interpole les vents, la temperature et le N**2 sur des isentropes
c       (en fait sur des iso-theta) lineairement en log(theta) =>
c       ucovteta, vcovteta, N2teta
c    3) on calcule la vorticite absolue sur des iso-theta => vorateta
c    4) on calcule la densite rho sur des iso-theta => rhoteta 
c
c       rhoteta = (T/theta)**(cp/R)*p0/(R*T)
c
c    5) on calcule la vorticite potentielle sur des iso-theta => PVteta
c
c       PVteta = (vorateta * N2 * theta)/(g * rhoteta) ! en PVU
c
c       NB: 1PVU=10**(-6) K*m**2/(s * kg)
c
c       PVteta =  vorateta * N2/(g**2 * rhoteta) ! en 1/(Pa*s)
c
c
c    *******************************************************************
c
c
c     Variables d'entree : ilon,ilev,pucov,pvcov,pteta,ztfi,zplay,zplev,nbteta,theta
c                       -> sur la grille dynamique
c     Variable de sortie : PVteta
c                       -> sur la grille physique 
c=======================================================================

c
c variables Input
c
      INTEGER ilon
      integer, intent(in):: ilev
      REAL, intent(in):: pvcov(iip1,jjm,ilev)
      REAL, intent(in):: pucov(iip1,jjp1,ilev)
      REAL, intent(in):: pteta(iip1,jjp1,ilev)
      REAL ztfi(ilon,ilev)
      REAL, intent(in):: zplay(ilon,ilev), zplev(ilon,ilev+1)
      INTEGER nbteta
      REAL theta(nbteta)
c
c variable Output
c
      REAL PVteta(ilon,nbteta)
c
c variables locales
c
      INTEGER i, j, l, ig0
      REAL SSUM
      REAL teta(ilon, ilev)
      REAL ptetau(ip1jmp1, ilev), ptetav(ip1jm, ilev)
      REAL ucovteta(ip1jmp1,ilev), vcovteta(ip1jm,ilev)
      REAL N2(ilon,ilev-1), N2teta(ilon,nbteta)
      REAL ztfiteta(ilon,nbteta)
      REAL rhoteta(ilon,nbteta)
      REAL vorateta(iip1,jjm,nbteta)
      REAL voratetafi(ilon,nbteta), vorpol(iim)
c
c
c projection teta sur la grille physique
c
      DO l=1,llm
       teta(1,l)   =  pteta(1,1,l)
       ig0         = 2
       DO j = 2, jjm
        DO i = 1, iim
         teta(ig0,l)    = pteta(i,j,l)
         ig0            = ig0 + 1
        ENDDO
       ENDDO
       teta(ig0,l)    = pteta(1,jjp1,l)
      ENDDO
c
c calcul pteta sur les grilles U et V
c
      DO l=1, llm
       DO j=1, jjp1
        DO i=1, iip1
         ig0=i+(j-1)*iip1
         ptetau(ig0,l)=pteta(i,j,l)
        ENDDO !i
       ENDDO !j
       DO j=1, jjm
        DO i=1, iip1
         ig0=i+(j-1)*iip1
         ptetav(ig0,l)=0.5*(pteta(i,j,l)+pteta(i,j+1,l))
        ENDDO !i
       ENDDO !j
      ENDDO !l
c
c projection pucov, pvcov sur une surface de theta constante
c
      DO l=1, nbteta
cIM 1rout CALL tetaleveli1j1(ip1jmp1,llm,.true.,ptetau,theta(l),
       CALL tetalevel(ip1jmp1,llm,.true.,ptetau,theta(l),
     .                pucov,ucovteta(:,l))
cIM 1rout CALL tetaleveli1j(ip1jm,llm,.true.,ptetav,theta(l),
       CALL tetalevel(ip1jm,llm,.true.,ptetav,theta(l),
     .                pvcov,vcovteta(:,l))
      ENDDO !l
c
c calcul vorticite absolue sur une iso-theta : vorateta
c
      CALL tourabs(nbteta,vcovteta,ucovteta,vorateta)
c
c projection vorateta sur la grille physique => voratetafi
c
      DO l=1,nbteta
       DO j=2,jjm
        ig0=1+(j-2)*iim
        DO i=1,iim
           voratetafi(ig0+i+1,l)
     &          = vorateta( i ,j-1,l) * alpha4_2d(i+1,j) +
     &          vorateta(i+1,j-1,l) * alpha1_2d(i+1,j) +
     &          vorateta(i  ,j  ,l) * alpha3_2d(i+1,j) +
     &          vorateta(i+1,j  ,l) * alpha2_2d(i+1,j)
        ENDDO
        voratetafi(ig0 +1,l) = voratetafi(ig0 +1+ iim,l)
       ENDDO
      ENDDO
c
      DO l=1,nbteta
       DO i=1,iim
        vorpol(i)  = vorateta(i,1,l)*aire_2d(i,1)
       ENDDO
       voratetafi(1,l)= SSUM(iim,vorpol,1)/apoln
      ENDDO
c
      DO l=1,nbteta
       DO i=1,iim
        vorpol(i)  = vorateta(i,jjm,l)* aire_2d(i,jjm +1)
       ENDDO
       voratetafi(ilon,l)= SSUM(iim,vorpol,1)/apols
      ENDDO
c 
c calcul N**2 sur la grille physique => N2
c
      DO l=1, llm-1 
       DO i=1, ilon
        N2(i,l) = (g**2 * zplay(i,l) * 
     &            (ztfi(i,l+1)-ztfi(i,l)) )/
     &            (R*ztfi(i,l)*ztfi(i,l)*
     &            (zplev(i,l)-zplev(i,l+1)) )+
     &            (g**2)/(ztfi(i,l)*CPP)
       ENDDO !i
      ENDDO !l
c
c calcul N2 sur une iso-theta => N2teta 
c
      DO l=1, nbteta
       CALL tetalevel(ilon,llm-1,.true.,teta,theta(l),
     &                N2,N2teta(:,l))
       CALL tetalevel(ilon,llm,.true.,teta,theta(l),
     &                ztfi,ztfiteta(:,l))
      ENDDO !l=1, nbteta
c
c calcul rho et PV sur une iso-theta : rhoteta, PVteta
c
      DO l=1, nbteta
       DO i=1, ilon
        rhoteta(i,l)=(ztfiteta(i,l)/theta(l))**(CPP/R)*
     &  (preff/(R*ztfiteta(i,l)))
c
c PVteta en PVU
c
        PVteta(i,l)=(theta(l)*g*voratetafi(i,l)*N2teta(i,l))/
     &              (g**2*rhoteta(i,l))
c
c PVteta en 1/(Pa*s)
c
        PVteta(i,l)=(voratetafi(i,l)*N2teta(i,l))/
     &              (g**2*rhoteta(i,l))
       ENDDO !i
      ENDDO !l
c
      RETURN
      END 
1	guez	3	SUBROUTINE PVtheta(ilon,ilev,pucov,pvcov,pteta,
2			& ztfi,zplay,zplev,
3			& nbteta,theta,PVteta)
4			use dimens_m
5			use paramet_m
6			use comconst
7	guez	66	use disvert_m
8	guez	3	use comgeom
9	guez	79	use tourabs_m, only: tourabs
10	guez	3	IMPLICIT none
11
12			c=======================================================================
13			c
14			c Auteur: I. Musat
15			c -------
16			c
17			c Objet:
18			c ------
19			c
20			c *******************************************************************
21			c Calcul de la vorticite potentielle PVteta sur des iso-theta selon
22			c la methodologie du NCEP/NCAR :
23			c 1) on calcule la stabilite statique N*2=g/T(dT/dz+g/cp) sur les
24			c niveaux du modele => N2
25			c 2) on interpole les vents, la temperature et le N**2 sur des isentropes
26			c (en fait sur des iso-theta) lineairement en log(theta) =>
27			c ucovteta, vcovteta, N2teta
28			c 3) on calcule la vorticite absolue sur des iso-theta => vorateta
29			c 4) on calcule la densite rho sur des iso-theta => rhoteta
30			c
31			c rhoteta = (T/theta)*(cp/R)p0/(R*T)
32			c
33			c 5) on calcule la vorticite potentielle sur des iso-theta => PVteta
34			c
35			c PVteta = (vorateta * N2 * theta)/(g * rhoteta) ! en PVU
36			c
37			c NB: 1PVU=10*(-6) Km*2/(s kg)
38			c
39			c PVteta = vorateta * N2/(g*2 rhoteta) ! en 1/(Pa*s)
40			c
41			c
42			c *******************************************************************
43			c
44			c
45			c Variables d'entree : ilon,ilev,pucov,pvcov,pteta,ztfi,zplay,zplev,nbteta,theta
46			c -> sur la grille dynamique
47			c Variable de sortie : PVteta
48			c -> sur la grille physique
49			c=======================================================================
50
51			c
52			c variables Input
53			c
54			INTEGER ilon
55			integer, intent(in):: ilev
56	guez	71	REAL, intent(in):: pvcov(iip1,jjm,ilev)
57	guez	70	REAL, intent(in):: pucov(iip1,jjp1,ilev)
58	guez	44	REAL, intent(in):: pteta(iip1,jjp1,ilev)
59	guez	3	REAL ztfi(ilon,ilev)
60	guez	10	REAL, intent(in):: zplay(ilon,ilev), zplev(ilon,ilev+1)
61	guez	3	INTEGER nbteta
62			REAL theta(nbteta)
63			c
64			c variable Output
65			c
66			REAL PVteta(ilon,nbteta)
67			c
68			c variables locales
69			c
70			INTEGER i, j, l, ig0
71			REAL SSUM
72			REAL teta(ilon, ilev)
73			REAL ptetau(ip1jmp1, ilev), ptetav(ip1jm, ilev)
74			REAL ucovteta(ip1jmp1,ilev), vcovteta(ip1jm,ilev)
75			REAL N2(ilon,ilev-1), N2teta(ilon,nbteta)
76			REAL ztfiteta(ilon,nbteta)
77			REAL rhoteta(ilon,nbteta)
78			REAL vorateta(iip1,jjm,nbteta)
79			REAL voratetafi(ilon,nbteta), vorpol(iim)
80			c
81			c
82			c projection teta sur la grille physique
83			c
84			DO l=1,llm
85			teta(1,l) = pteta(1,1,l)
86			ig0 = 2
87			DO j = 2, jjm
88			DO i = 1, iim
89			teta(ig0,l) = pteta(i,j,l)
90			ig0 = ig0 + 1
91			ENDDO
92			ENDDO
93			teta(ig0,l) = pteta(1,jjp1,l)
94			ENDDO
95			c
96			c calcul pteta sur les grilles U et V
97			c
98			DO l=1, llm
99			DO j=1, jjp1
100			DO i=1, iip1
101			ig0=i+(j-1)*iip1
102			ptetau(ig0,l)=pteta(i,j,l)
103			ENDDO !i
104			ENDDO !j
105			DO j=1, jjm
106			DO i=1, iip1
107			ig0=i+(j-1)*iip1
108			ptetav(ig0,l)=0.5*(pteta(i,j,l)+pteta(i,j+1,l))
109			ENDDO !i
110			ENDDO !j
111			ENDDO !l
112			c
113			c projection pucov, pvcov sur une surface de theta constante
114			c
115			DO l=1, nbteta
116			cIM 1rout CALL tetaleveli1j1(ip1jmp1,llm,.true.,ptetau,theta(l),
117			CALL tetalevel(ip1jmp1,llm,.true.,ptetau,theta(l),
118			. pucov,ucovteta(:,l))
119			cIM 1rout CALL tetaleveli1j(ip1jm,llm,.true.,ptetav,theta(l),
120			CALL tetalevel(ip1jm,llm,.true.,ptetav,theta(l),
121			. pvcov,vcovteta(:,l))
122			ENDDO !l
123			c
124			c calcul vorticite absolue sur une iso-theta : vorateta
125			c
126			CALL tourabs(nbteta,vcovteta,ucovteta,vorateta)
127			c
128			c projection vorateta sur la grille physique => voratetafi
129			c
130			DO l=1,nbteta
131			DO j=2,jjm
132			ig0=1+(j-2)*iim
133			DO i=1,iim
134			voratetafi(ig0+i+1,l)
135			& = vorateta( i ,j-1,l) * alpha4_2d(i+1,j) +
136			& vorateta(i+1,j-1,l) * alpha1_2d(i+1,j) +
137			& vorateta(i ,j ,l) * alpha3_2d(i+1,j) +
138			& vorateta(i+1,j ,l) * alpha2_2d(i+1,j)
139			ENDDO
140			voratetafi(ig0 +1,l) = voratetafi(ig0 +1+ iim,l)
141			ENDDO
142			ENDDO
143			c
144			DO l=1,nbteta
145			DO i=1,iim
146			vorpol(i) = vorateta(i,1,l)*aire_2d(i,1)
147			ENDDO
148			voratetafi(1,l)= SSUM(iim,vorpol,1)/apoln
149			ENDDO
150			c
151			DO l=1,nbteta
152			DO i=1,iim
153			vorpol(i) = vorateta(i,jjm,l)* aire_2d(i,jjm +1)
154			ENDDO
155			voratetafi(ilon,l)= SSUM(iim,vorpol,1)/apols
156			ENDDO
157			c
158			c calcul N**2 sur la grille physique => N2
159			c
160			DO l=1, llm-1
161			DO i=1, ilon
162			N2(i,l) = (g*2 zplay(i,l) *
163			& (ztfi(i,l+1)-ztfi(i,l)) )/
164			& (Rztfi(i,l)ztfi(i,l)*
165			& (zplev(i,l)-zplev(i,l+1)) )+
166			& (g*2)/(ztfi(i,l)CPP)
167			ENDDO !i
168			ENDDO !l
169			c
170			c calcul N2 sur une iso-theta => N2teta
171			c
172			DO l=1, nbteta
173			CALL tetalevel(ilon,llm-1,.true.,teta,theta(l),
174			& N2,N2teta(:,l))
175			CALL tetalevel(ilon,llm,.true.,teta,theta(l),
176			& ztfi,ztfiteta(:,l))
177			ENDDO !l=1, nbteta
178			c
179			c calcul rho et PV sur une iso-theta : rhoteta, PVteta
180			c
181			DO l=1, nbteta
182			DO i=1, ilon
183			rhoteta(i,l)=(ztfiteta(i,l)/theta(l))*(CPP/R)
184			& (preff/(R*ztfiteta(i,l)))
185			c
186			c PVteta en PVU
187			c
188			PVteta(i,l)=(theta(l)gvoratetafi(i,l)*N2teta(i,l))/
189			& (g*2rhoteta(i,l))
190			c
191			c PVteta en 1/(Pa*s)
192			c
193			PVteta(i,l)=(voratetafi(i,l)*N2teta(i,l))/
194			& (g*2rhoteta(i,l))
195			ENDDO !i
196			ENDDO !l
197			c
198			RETURN
199			END