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

Diff of /trunk/dyn3d/PVtheta.f

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

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