source: LMDZ5/branches/LMDZ5_AR5/libf/dyn3d/PVtheta.F @ 4116

Last change on this file since 4116 was 649, checked in by lmdzadmin, 19 years ago

Version erronee de la routine originelle
IM/LF

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