trunk/dyn3d/PVtheta.f

module pvtheta_m

  IMPLICIT NONE

contains

  SUBROUTINE pvtheta(ilon, ilev, pucov, pvcov, pteta, ztfi, zplay, zplev, &
       nbteta, theta, pvteta)

    USE comconst, ONLY: cpp, g, r
    USE comgeom, ONLY: aire_2d, alpha1_2d, alpha2_2d, alpha3_2d, alpha4_2d, &
         apoln, apols
    USE dimens_m, ONLY: iim, jjm, llm
    USE disvert_m, ONLY: preff
    USE paramet_m, ONLY: iip1, ip1jm, ip1jmp1, jjp1
    USE tourabs_m, ONLY: tourabs

    ! Author: I. Musat

    ! Calcul de la vorticite potentielle PVteta sur des iso-theta selon
    ! la methodologie du NCEP/NCAR :

    ! 1) on calcule n2

    ! 2) on interpole les vents, la temperature et le N2 sur des
    ! isentropes (en fait sur des iso-theta) lineairement en
    ! log(theta), d'où ucovteta, vcovteta, N2teta

    ! 3) on calcule vorateta
    ! 4) on calcule rhoteta
    ! 5) on calcule PVteta

    INTEGER, INTENT(IN):: ilon
    INTEGER, INTENT(IN):: ilev

    ! pucov, pvcov, pteta, ztfi, zplay, zplev sur la grille dynamique
    REAL, INTENT(IN):: pucov(iip1, jjp1, ilev)
    REAL, INTENT(IN):: pvcov(iip1, jjm, ilev)
    REAL, INTENT(IN):: pteta(iip1, jjp1, ilev)
    REAL, INTENT(IN):: ztfi(ilon, ilev)
    REAL, INTENT(IN):: zplay(ilon, ilev), zplev(ilon, ilev + 1)

    INTEGER, INTENT(IN):: nbteta
    REAL, INTENT(IN):: theta(nbteta) ! sur la grille dynamique

    REAL, INTENT(out):: pvteta(ilon, nbteta) 
    ! vorticité potentielle sur des iso-theta, sur la grille physique,
    ! en Pa-1 s-1
    ! PVteta = vorateta * N2 / (g**2 * rhoteta)

    ! Local:

    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) ! stabilité statique sur les niveaux du modèle
    ! N**2 = g / T * (dT / dz + g / cp)

    real n2teta(ilon, nbteta) ! N**2 sur une iso-theta
    REAL ztfiteta(ilon, nbteta)

    REAL rhoteta(ilon, nbteta) ! densite sur des iso-theta
    ! rhoteta = (T / theta)**(cp / R) * p0 / (R * T)

    REAL vorateta(iip1, jjm, nbteta) ! vorticite absolue sur des iso-theta
    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
       END DO
       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
       END DO
    END DO

    ! projection pucov, pvcov sur une surface de theta constante
    DO l = 1, nbteta
       CALL tetalevel(ip1jmp1, llm, .TRUE., ptetau, theta(l), pucov, &
            ucovteta(:, l))
       CALL tetalevel(ip1jm, llm, .TRUE., ptetav, theta(l), pvcov, &
            vcovteta(:, l))
    END DO

    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

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

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

    DO l = 1, nbteta
       DO i = 1, ilon
          rhoteta(i, l) = (ztfiteta(i, l) / theta(l))**(cpp / r) &
               * (preff / (r * ztfiteta(i, l)))
          pvteta(i, l) = (voratetafi(i, l) * n2teta(i, l)) &
               / (g**2 * rhoteta(i, l))
       END DO
    END DO

  END SUBROUTINE pvtheta

end module pvtheta_m
1	guez	88	module pvtheta_m
2
3	guez	81	IMPLICIT NONE
4	guez	3
5	guez	88	contains
6	guez	3
7	guez	88	SUBROUTINE pvtheta(ilon, ilev, pucov, pvcov, pteta, ztfi, zplay, zplev, &
8			nbteta, theta, pvteta)
9	guez	81
10	guez	88	USE comconst, ONLY: cpp, g, r
11			USE comgeom, ONLY: aire_2d, alpha1_2d, alpha2_2d, alpha3_2d, alpha4_2d, &
12			apoln, apols
13			USE dimens_m, ONLY: iim, jjm, llm
14			USE disvert_m, ONLY: preff
15			USE paramet_m, ONLY: iip1, ip1jm, ip1jmp1, jjp1
16			USE tourabs_m, ONLY: tourabs
17	guez	81
18	guez	88	! Author: I. Musat
19	guez	81
20	guez	88	! Calcul de la vorticite potentielle PVteta sur des iso-theta selon
21			! la methodologie du NCEP/NCAR :
22	guez	81
23	guez	88	! 1) on calcule n2
24	guez	81
25	guez	88	! 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	guez	81
29	guez	88	! 3) on calcule vorateta
30			! 4) on calcule rhoteta
31			! 5) on calcule PVteta
32	guez	81
33	guez	88	INTEGER, INTENT(IN):: ilon
34			INTEGER, INTENT(IN):: ilev
35	guez	81
36	guez	88	! 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	guez	81
43	guez	88	INTEGER, INTENT(IN):: nbteta
44			REAL, INTENT(IN):: theta(nbteta) ! sur la grille dynamique
45	guez	81
46	guez	88	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	guez	81
51	guez	88	! Local:
52	guez	81
53	guez	88	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	guez	81
59	guez	88	REAL n2(ilon, ilev - 1) ! stabilité statique sur les niveaux du modèle
60			! N*2 = g / T (dT / dz + g / cp)
61	guez	81
62	guez	88	real n2teta(ilon, nbteta) ! N**2 sur une iso-theta
63			REAL ztfiteta(ilon, nbteta)
64	guez	81
65	guez	88	REAL rhoteta(ilon, nbteta) ! densite sur des iso-theta
66			! rhoteta = (T / theta)*(cp / R) p0 / (R * T)
67	guez	81
68	guez	88	REAL vorateta(iip1, jjm, nbteta) ! vorticite absolue sur des iso-theta
69			REAL voratetafi(ilon, nbteta), vorpol(iim)
70	guez	81
71	guez	88	!-------------------------------------------------------------------------
72	guez	81
73	guez	88	! 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
85	guez	81
86	guez	88	! calcul pteta sur les grilles U et V
87			DO l = 1, llm
88			DO j = 1, jjp1
89			DO i = 1, iip1
90			ig0 = i + (j - 1) * iip1
91			ptetau(ig0, l) = pteta(i, j, l)
92			END DO
93			END DO
94			DO j = 1, jjm
95			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	guez	81
102	guez	88	! 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	guez	81	END DO
109
110	guez	88	CALL tourabs(nbteta, vcovteta, ucovteta, vorateta)
111	guez	81
112	guez	88	! projection vorateta sur la grille physique => voratetafi
113	guez	81
114	guez	88	DO l = 1, nbteta
115			DO j = 2, jjm
116			ig0 = 1 + (j - 2) * iim
117			DO i = 1, iim
118			voratetafi(ig0 + i + 1, l) = vorateta(i, j - 1, l) &
119			* alpha4_2d(i + 1, j) + vorateta(i + 1, j - 1, l) &
120			* alpha1_2d(i + 1, j) + vorateta(i, j, l) &
121			* alpha3_2d(i + 1, j) + vorateta(i + 1, j, l) &
122			* alpha2_2d(i + 1, j)
123			END DO
124			voratetafi(ig0 + 1, l) = voratetafi(ig0 + 1 + iim, l)
125			END DO
126			END DO
127	guez	81
128	guez	88	DO l = 1, nbteta
129			DO i = 1, iim
130			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
134	guez	81
135	guez	88	DO l = 1, nbteta
136			DO i = 1, iim
137			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
141	guez	81
142	guez	88	DO l = 1, llm - 1
143			DO i = 1, ilon
144			n2(i, l) = (g*2 zplay(i, l) * (ztfi(i, l + 1) - ztfi(i, l))) &
145			/ (r * ztfi(i, l) * ztfi(i, l) &
146			* (zplev(i, l) - zplev(i, l + 1))) + (g*2) / (ztfi(i, l) cpp)
147			END DO
148	guez	81	END DO
149
150	guez	88	! calcul N2teta
151			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	guez	81	END DO
155
156	guez	88	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	guez	81	END DO
164
165	guez	88	END SUBROUTINE pvtheta
166	guez	81
167	guez	88	end module pvtheta_m