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