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

  ! =======================================================================

  ! 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

  REAL pvteta(ilon, nbteta)

  ! 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)


  ! projection teta sur la grille physique

  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
      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

  ! projection pucov, pvcov sur une surface de theta constante

  DO l = 1, nbteta
    ! IM 1rout CALL tetaleveli1j1(ip1jmp1,llm,.true.,ptetau,theta(l),
    CALL tetalevel(ip1jmp1, llm, .TRUE., ptetau, theta(l), pucov, &
      ucovteta(:,l))
    ! IM 1rout CALL tetaleveli1j(ip1jm,llm,.true.,ptetav,theta(l),
    CALL tetalevel(ip1jm, llm, .TRUE., ptetav, theta(l), pvcov, &
      vcovteta(:,l))
  END DO !l

  ! calcul vorticite absolue sur une iso-theta : vorateta

  CALL tourabs(nbteta, vcovteta, ucovteta, vorateta)

  ! projection vorateta sur la grille physique => voratetafi

  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)
      END DO
      voratetafi(ig0+1, l) = voratetafi(ig0+1+iim, l)
    END DO
  END DO

  DO l = 1, nbteta
    DO i = 1, iim
      vorpol(i) = vorateta(i, 1, l)*aire_2d(i, 1)
    END DO
    voratetafi(1, l) = ssum(iim, vorpol, 1)/apoln
  END DO

  DO l = 1, nbteta
    DO i = 1, iim
      vorpol(i) = vorateta(i, jjm, l)*aire_2d(i, jjm+1)
    END DO
    voratetafi(ilon, l) = ssum(iim, vorpol, 1)/apols
  END DO

  ! calcul N**2 sur la grille physique => N2

  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)
    END DO !i
  END DO !l

  ! 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

      pvteta(i, l) = (theta(l)*g*voratetafi(i,l)*n2teta(i,l))/ &
        (g**2*rhoteta(i,l))

      ! PVteta en 1/(Pa*s)

      pvteta(i, l) = (voratetafi(i,l)*n2teta(i,l))/(g**2*rhoteta(i,l))
    END DO !i
  END DO !l

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