/[lmdze]/trunk/dyn3d/PVtheta.f90
ViewVC logotype

Diff of /trunk/dyn3d/PVtheta.f90

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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