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	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
11	! =======================================================================
12
13	! 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