16 |
! tangente hyperbolique |
! tangente hyperbolique |
17 |
! calcul des coefficients (cu_2d, cv_2d, 1./cu_2d**2, 1./cv_2d**2) |
! calcul des coefficients (cu_2d, cv_2d, 1./cu_2d**2, 1./cv_2d**2) |
18 |
|
|
19 |
! les coef. ( cu_2d, cv_2d ) permettent de passer des vitesses naturelles |
! les coef. (cu_2d, cv_2d) permettent de passer des vitesses naturelles |
20 |
! aux vitesses covariantes et contravariantes , ou vice-versa ... |
! aux vitesses covariantes et contravariantes, ou vice-versa |
21 |
|
|
22 |
! on a : u (covariant) = cu_2d * u (naturel) , u(contrav)= u(nat)/cu_2d |
! on a : |
23 |
! v (covariant) = cv_2d * v (naturel) , v(contrav)= v(nat)/cv_2d |
! u (covariant) = cu_2d * u (naturel), u(contrav)= u(nat)/cu_2d |
24 |
|
! v (covariant) = cv_2d * v (naturel), v(contrav)= v(nat)/cv_2d |
25 |
! on en tire : u(covariant) = cu_2d * cu_2d * u(contravariant) |
|
26 |
! v(covariant) = cv_2d * cv_2d * v(contravariant) |
! on en tire : |
27 |
|
! u(covariant) = cu_2d * cu_2d * u(contravariant) |
28 |
! on a l'application ( x(X) , y(Y) ) avec - im/2 +1 < X < im/2 |
! v(covariant) = cv_2d * cv_2d * v(contravariant) |
29 |
! = = |
|
30 |
! et - jm/2 < Y < jm/2 |
! on a l'application (x(X), y(Y)) avec - im/2 +1 <= X <= im/2 |
31 |
! = = |
! et - jm/2 <= Y <= jm/2 |
32 |
|
|
33 |
! . x est la longitude du point en radians . |
! x est la longitude du point en radians. |
34 |
! . y est la latitude du point en radians . |
! y est la latitude du point en radians. |
35 |
! . . |
! |
36 |
! . on a : cu_2d(i, j) = rad * COS(y) * dx/dX . |
! on a : cu_2d(i, j) = rad * cos(y) * dx/dX |
37 |
! . cv( j ) = rad * dy/dY . |
! cv(j) = rad * dy/dY |
38 |
! . aire_2d(i, j) = cu_2d(i, j) * cv(j) . |
! aire_2d(i, j) = cu_2d(i, j) * cv(j) |
39 |
! . . |
! |
40 |
! . y, dx/dX, dy/dY calcules aux points concernes . |
! y, dx/dX, dy/dY calcules aux points concernes |
41 |
! , |
! cv, bien que dependant de j uniquement, sera ici indice aussi en i |
42 |
! cv , bien que dependant de j uniquement, sera ici indice aussi en i |
! pour un adressage plus facile en ij. |
43 |
! pour un adressage plus facile en ij . |
|
44 |
|
! aux points u et v, |
45 |
! ************** aux points u et v , ***************** |
! xprimu et xprimv sont respectivement les valeurs de dx/dX |
46 |
! xprimu et xprimv sont respectivement les valeurs de dx/dX |
! yprimu et yprimv sont respectivement les valeurs de dy/dY |
47 |
! yprimu et yprimv . . . . . . . . . . . dy/dY |
! rlatu et rlatv sont respectivement les valeurs de la latitude |
48 |
! rlatu et rlatv . . . . . . . . . . .la latitude |
! cvu et cv_2d sont respectivement les valeurs de cv_2d |
49 |
! cvu et cv_2d . . . . . . . . . . . cv_2d |
|
50 |
|
! aux points u, v, scalaires, et z |
51 |
! ************** aux points u, v, scalaires, et z **************** |
! cu_2d, cuv, cuscal, cuz sont respectivement les valeurs de cu_2d |
52 |
! cu_2d, cuv, cuscal, cuz sont respectiv. les valeurs de cu_2d |
! Cf. "inigeom.txt". |
|
|
|
|
! Exemple de distribution de variables sur la grille dans le |
|
|
! domaine de travail ( X, Y ) . |
|
|
! DX=DY= 1 |
|
|
|
|
|
! + represente un point scalaire ( p.exp la pression ) |
|
|
! > represente la composante zonale du vent |
|
|
! V represente la composante meridienne du vent |
|
|
! o represente la vorticite |
|
|
|
|
|
! ---- , car aux poles , les comp.zonales covariantes sont nulles |
|
|
|
|
|
! i -> |
|
|
|
|
|
! 1 2 3 4 5 6 7 8 |
|
|
! j |
|
|
! v 1 + ---- + ---- + ---- + ---- + ---- + ---- + ---- + -- |
|
|
|
|
|
! V o V o V o V o V o V o V o V o |
|
|
|
|
|
! 2 + > + > + > + > + > + > + > + > |
|
|
|
|
|
! V o V o V o V o V o V o V o V o |
|
|
|
|
|
! 3 + > + > + > + > + > + > + > + > |
|
|
|
|
|
! V o V o V o V o V o V o V o V o |
|
|
|
|
|
! 4 + > + > + > + > + > + > + > + > |
|
|
|
|
|
! V o V o V o V o V o V o V o V o |
|
|
|
|
|
! 5 + ---- + ---- + ---- + ---- + ---- + ---- + ---- + -- |
|
|
|
|
|
! Ci-dessus, on voit que le nombre de pts.en longitude est egal |
|
|
! a IM = 8 |
|
|
! De meme , le nombre d'intervalles entre les 2 poles est egal |
|
|
! a JM = 4 |
|
|
|
|
|
! Les points scalaires ( + ) correspondent donc a des valeurs |
|
|
! entieres de i ( 1 a IM ) et de j ( 1 a JM +1 ) . |
|
|
|
|
|
! Les vents U ( > ) correspondent a des valeurs semi- |
|
|
! entieres de i ( 1+ 0.5 a IM+ 0.5) et entieres de j ( 1 a JM+1) |
|
|
|
|
|
! Les vents V ( V ) correspondent a des valeurs entieres |
|
|
! de i ( 1 a IM ) et semi-entieres de j ( 1 +0.5 a JM +0.5) |
|
53 |
|
|
54 |
USE dimens_m, ONLY : iim, jjm |
USE dimens_m, ONLY : iim, jjm |
55 |
USE paramet_m, ONLY : iip1, jjp1 |
USE paramet_m, ONLY : iip1, jjp1 |
92 |
real aireij2_2d(iim + 1, jjm + 1) |
real aireij2_2d(iim + 1, jjm + 1) |
93 |
real aireij3_2d(iim + 1, jjm + 1), aireij4_2d(iim + 1, jjm + 1) |
real aireij3_2d(iim + 1, jjm + 1), aireij4_2d(iim + 1, jjm + 1) |
94 |
real airuscv2_2d(iim + 1, jjm) |
real airuscv2_2d(iim + 1, jjm) |
95 |
real airvscu2_2d(iim + 1, jjm), aiuscv2gam_2d(iim + 1, jjm) |
real airvscu2_2d(iim + 1, jjm), aiuscv2gam_2d(iim + 1, jjm) |
96 |
real aivscu2gam_2d(iim + 1, jjm) |
real aivscu2gam_2d(iim + 1, jjm) |
97 |
|
|
98 |
!------------------------------------------------------------------ |
!------------------------------------------------------------------ |
121 |
|
|
122 |
WRITE (6, 990) |
WRITE (6, 990) |
123 |
|
|
124 |
IF ( .NOT. fxyhypb) THEN |
IF (.NOT. fxyhypb) THEN |
125 |
IF (ysinus) THEN |
IF (ysinus) THEN |
126 |
print *, ' *** Inigeom , Y = Sinus ( Latitude ) *** ' |
print *, ' Inigeom, Y = Sinus (Latitude) ' |
127 |
|
! utilisation de f(x, y) avec y = sinus de la latitude |
|
! utilisation de f(x, y ) avec y = sinus de la latitude ... |
|
|
|
|
128 |
CALL fxysinus(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, & |
CALL fxysinus(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, & |
129 |
rlatu2, yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, & |
rlatu2, yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, & |
130 |
xprimm025, rlonp025, xprimp025) |
xprimm025, rlonp025, xprimp025) |
131 |
ELSE |
ELSE |
132 |
print *, '*** Inigeom , Y = Latitude , der. sinusoid . ***' |
print *, 'Inigeom, Y = Latitude, der. sinusoid .' |
133 |
! utilisation de f(x, y) a tangente sinusoidale , y etant la latit |
! utilisation de f(x, y) a tangente sinusoidale, y etant la latit |
134 |
|
|
135 |
pxo = clon*pi/180. |
pxo = clon*pi/180. |
136 |
pyo = 2.*clat*pi/180. |
pyo = 2.*clat*pi/180. |
137 |
|
|
138 |
! determination de transx ( pour le zoom ) par Newton-Raphson . |
! determination de transx (pour le zoom) par Newton-Raphson |
139 |
|
|
140 |
itmax = 10 |
itmax = 10 |
141 |
eps = .1E-7 |
eps = .1E-7 |
174 |
rlonp025, xprimp025) |
rlonp025, xprimp025) |
175 |
END IF |
END IF |
176 |
ELSE |
ELSE |
177 |
! .... Utilisation de fxyhyper , f(x, y) a derivee tangente hyperbol. |
! Utilisation de fxyhyper, f(x, y) à dérivée tangente hyperbolique |
178 |
print *, '*** Inigeom , Y = Latitude , der.tg. hyperbolique ***' |
print *, 'Inigeom, Y = Latitude, dérivée tangente hyperbolique' |
179 |
CALL fxyhyper(clat, grossismy, dzoomy, tauy, clon, grossismx, dzoomx, & |
CALL fxyhyper(clat, grossismy, dzoomy, tauy, clon, grossismx, dzoomx, & |
180 |
taux, rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, & |
taux, rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, & |
181 |
yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, & |
yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, & |
185 |
rlatu(1) = asin(1.) |
rlatu(1) = asin(1.) |
186 |
rlatu(jjp1) = -rlatu(1) |
rlatu(jjp1) = -rlatu(1) |
187 |
|
|
188 |
! .... calcul aux poles .... |
! calcul aux poles |
189 |
|
|
190 |
yprimu(1) = 0. |
yprimu(1) = 0. |
191 |
yprimu(jjp1) = 0. |
yprimu(jjp1) = 0. |
192 |
|
|
193 |
un4rad2 = 0.25*rad*rad |
un4rad2 = 0.25*rad*rad |
194 |
|
|
195 |
! calcul des aires ( aire_2d, aireu_2d, airev_2d, 1./aire_2d, 1./airez ) |
! calcul des aires (aire_2d, aireu_2d, airev_2d, 1./aire_2d, 1./airez) |
196 |
! - et de fext_2d , force de coriolis extensive . |
! - et de fext_2d, force de coriolis extensive |
197 |
|
|
198 |
! A 1 point scalaire P (i, j) de la grille, reguliere en (X, Y) , sont |
! A 1 point scalaire P (i, j) de la grille, reguliere en (X, Y), sont |
199 |
! affectees 4 aires entourant P , calculees respectivement aux points |
! affectees 4 aires entourant P, calculees respectivement aux points |
200 |
! ( i + 1/4, j - 1/4 ) : aireij1_2d (i, j) |
! (i + 1/4, j - 1/4) : aireij1_2d (i, j) |
201 |
! ( i + 1/4, j + 1/4 ) : aireij2_2d (i, j) |
! (i + 1/4, j + 1/4) : aireij2_2d (i, j) |
202 |
! ( i - 1/4, j + 1/4 ) : aireij3_2d (i, j) |
! (i - 1/4, j + 1/4) : aireij3_2d (i, j) |
203 |
! ( i - 1/4, j - 1/4 ) : aireij4_2d (i, j) |
! (i - 1/4, j - 1/4) : aireij4_2d (i, j) |
204 |
|
|
205 |
! , |
!, |
206 |
! Les cotes de chacun de ces 4 carres etant egaux a 1/2 suivant (X, Y). |
! Les cotes de chacun de ces 4 carres etant egaux a 1/2 suivant (X, Y). |
207 |
! Chaque aire centree en 1 point scalaire P(i, j) est egale a la somme |
! Chaque aire centree en 1 point scalaire P(i, j) est egale a la somme |
208 |
! des 4 aires aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d qui sont |
! des 4 aires aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d qui sont |
209 |
! affectees au |
! affectees au |
210 |
! point (i, j) . |
! point (i, j). |
211 |
! On definit en outre les coefficients alpha comme etant egaux a |
! On definit en outre les coefficients alpha comme etant egaux a |
212 |
! (aireij / aire_2d), c.a.d par exp. |
! (aireij / aire_2d), c.a.d par exp. |
213 |
! alpha1_2d(i, j)=aireij1_2d(i, j)/aire_2d(i, j) |
! alpha1_2d(i, j)=aireij1_2d(i, j)/aire_2d(i, j) |
214 |
|
|
215 |
! De meme, toute aire centree en 1 point U est egale a la somme des |
! De meme, toute aire centree en 1 point U est egale a la somme des |
216 |
! 4 aires aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d entourant |
! 4 aires aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d entourant |
217 |
! le point U. |
! le point U. |
218 |
! Idem pour airev_2d, airez . |
! Idem pour airev_2d, airez. |
219 |
|
|
220 |
! On a , pour chaque maille : dX = dY = 1 |
! On a, pour chaque maille : dX = dY = 1 |
221 |
|
|
222 |
! . V |
! V |
223 |
|
|
224 |
! aireij4_2d . . aireij1_2d |
! aireij4_2d . . aireij1_2d |
225 |
|
|
226 |
! U . . P . U |
! U . . P . U |
227 |
|
|
228 |
! aireij3_2d . . aireij2_2d |
! aireij3_2d . . aireij2_2d |
229 |
|
|
230 |
! . V |
! V |
231 |
|
|
232 |
! Calcul des 4 aires elementaires aireij1_2d, aireij2_2d, |
! Calcul des 4 aires elementaires aireij1_2d, aireij2_2d, |
233 |
! aireij3_2d, aireij4_2d |
! aireij3_2d, aireij4_2d |
234 |
! qui entourent chaque aire_2d(i, j) , ainsi que les 4 elongations |
! qui entourent chaque aire_2d(i, j), ainsi que les 4 elongations |
235 |
! elementaires |
! elementaires |
236 |
! cuij et les 4 elongat. cvij qui sont calculees aux memes |
! cuij et les 4 elongat. cvij qui sont calculees aux memes |
237 |
! endroits que les aireij . |
! endroits que les aireij. |
238 |
|
|
239 |
! ....... do 35 : boucle sur les jjm + 1 latitudes ..... |
! do 35 : boucle sur les jjm + 1 latitudes |
240 |
|
|
241 |
DO j = 1, jjp1 |
DO j = 1, jjp1 |
242 |
|
|
333 |
|
|
334 |
END IF |
END IF |
335 |
|
|
336 |
! ........ periodicite ............ |
! periodicite |
337 |
|
|
338 |
cvij1(iip1, j) = cvij1(1, j) |
cvij1(iip1, j) = cvij1(1, j) |
339 |
cvij2(iip1, j) = cvij2(1, j) |
cvij2(iip1, j) = cvij2(1, j) |
411 |
|
|
412 |
END DO |
END DO |
413 |
|
|
414 |
! ..... Calcul des elongations cu_2d, cv_2d, cvu ......... |
! Calcul des elongations cu_2d, cv_2d, cvu |
415 |
|
|
416 |
DO j = 1, jjm |
DO j = 1, jjm |
417 |
DO i = 1, iim |
DO i = 1, iim |
459 |
cuscvugam_2d(iip1, j) = cuscvugam_2d(1, j) |
cuscvugam_2d(iip1, j) = cuscvugam_2d(1, j) |
460 |
END DO |
END DO |
461 |
|
|
462 |
! .... calcul aux poles .... |
! calcul aux poles |
463 |
|
|
464 |
DO i = 1, iip1 |
DO i = 1, iip1 |
465 |
cu_2d(i, 1) = 0. |
cu_2d(i, 1) = 0. |
489 |
aiuscv2gam_2d(iip1, j) = aiuscv2gam_2d(1, j) |
aiuscv2gam_2d(iip1, j) = aiuscv2gam_2d(1, j) |
490 |
END DO |
END DO |
491 |
|
|
492 |
! calcul des aires aux poles : |
! calcul des aires aux poles : |
493 |
|
|
494 |
apoln = sum(aire_2d(:iim, 1)) |
apoln = sum(aire_2d(:iim, 1)) |
495 |
apols = sum(aire_2d(:iim, jjp1)) |
apols = sum(aire_2d(:iim, jjp1)) |
498 |
unsapolnga2 = 1./(apoln**(-gamdi_h)) |
unsapolnga2 = 1./(apoln**(-gamdi_h)) |
499 |
unsapolsga2 = 1./(apols**(-gamdi_h)) |
unsapolsga2 = 1./(apols**(-gamdi_h)) |
500 |
|
|
501 |
! changement F. Hourdin calcul conservatif pour fext_2d |
! changement F. Hourdin calcul conservatif pour fext_2d |
502 |
! constang_2d contient le produit a * cos ( latitude ) * omega |
! constang_2d contient le produit a * cos (latitude) * omega |
503 |
|
|
504 |
DO i = 1, iim |
DO i = 1, iim |
505 |
constang_2d(i, 1) = 0. |
constang_2d(i, 1) = 0. |
513 |
constang_2d(i, jjp1) = 0. |
constang_2d(i, jjp1) = 0. |
514 |
END DO |
END DO |
515 |
|
|
516 |
! periodicite en longitude |
! periodicite en longitude |
517 |
|
|
518 |
DO j = 1, jjm |
DO j = 1, jjm |
519 |
fext_2d(iip1, j) = fext_2d(1, j) |
fext_2d(iip1, j) = fext_2d(1, j) |
524 |
|
|
525 |
! fin du changement |
! fin du changement |
526 |
|
|
527 |
print *, ' *** Coordonnees de la grille *** ' |
print *, ' Coordonnees de la grille ' |
528 |
WRITE (6, 995) |
WRITE (6, 995) |
529 |
|
|
530 |
print *, ' LONGITUDES aux pts. V ( degres ) ' |
print *, ' LONGITUDES aux pts. V (degres) ' |
531 |
WRITE (6, 995) |
WRITE (6, 995) |
532 |
DO i = 1, iip1 |
DO i = 1, iip1 |
533 |
rlonvv(i) = rlonv(i)*180./pi |
rlonvv(i) = rlonv(i)*180./pi |
535 |
WRITE (6, 400) rlonvv |
WRITE (6, 400) rlonvv |
536 |
|
|
537 |
WRITE (6, 995) |
WRITE (6, 995) |
538 |
print *, ' LATITUDES aux pts. V ( degres ) ' |
print *, ' LATITUDES aux pts. V (degres) ' |
539 |
WRITE (6, 995) |
WRITE (6, 995) |
540 |
DO i = 1, jjm |
DO i = 1, jjm |
541 |
rlatuu(i) = rlatv(i)*180./pi |
rlatuu(i) = rlatv(i)*180./pi |
546 |
rlonvv(i) = rlonu(i)*180./pi |
rlonvv(i) = rlonu(i)*180./pi |
547 |
END DO |
END DO |
548 |
WRITE (6, 995) |
WRITE (6, 995) |
549 |
print *, ' LONGITUDES aux pts. U ( degres ) ' |
print *, ' LONGITUDES aux pts. U (degres) ' |
550 |
WRITE (6, 995) |
WRITE (6, 995) |
551 |
WRITE (6, 400) rlonvv |
WRITE (6, 400) rlonvv |
552 |
WRITE (6, 995) |
WRITE (6, 995) |
553 |
|
|
554 |
print *, ' LATITUDES aux pts. U ( degres ) ' |
print *, ' LATITUDES aux pts. U (degres) ' |
555 |
WRITE (6, 995) |
WRITE (6, 995) |
556 |
DO i = 1, jjp1 |
DO i = 1, jjp1 |
557 |
rlatuu(i) = rlatu(i)*180./pi |
rlatuu(i) = rlatu(i)*180./pi |