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 |
|
|
26 |
|
! on en tire : |
27 |
|
! u(covariant) = cu_2d * cu_2d * u(contravariant) |
28 |
|
! v(covariant) = cv_2d * cv_2d * v(contravariant) |
29 |
|
|
30 |
|
! 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. |
34 |
|
! y est la latitude du point en radians. |
35 |
|
! |
36 |
|
! on a : cu_2d(i, j) = rad * cos(y) * dx/dX |
37 |
|
! cv(j) = rad * dy/dY |
38 |
|
! aire_2d(i, j) = cu_2d(i, j) * cv(j) |
39 |
|
! |
40 |
|
! y, dx/dX, dy/dY calcules aux points concernes |
41 |
|
! cv, bien que dependant de j uniquement, sera ici indice aussi en i |
42 |
|
! pour un adressage plus facile en ij. |
43 |
|
|
44 |
|
! aux points u et v, |
45 |
|
! xprimu et xprimv sont respectivement les valeurs de dx/dX |
46 |
|
! yprimu et yprimv sont respectivement les valeurs de dy/dY |
47 |
|
! rlatu et rlatv sont respectivement les valeurs de la latitude |
48 |
|
! cvu et cv_2d sont respectivement les valeurs de cv_2d |
49 |
|
|
50 |
|
! aux points u, v, scalaires, et z |
51 |
|
! cu_2d, cuv, cuscal, cuz sont respectivement les valeurs de cu_2d |
52 |
|
! Cf. "inigeom.txt". |
53 |
|
|
54 |
! on en tire : u(covariant) = cu_2d * cu_2d * u(contravariant) |
USE comconst, ONLY : g, omeg, rad |
|
! v(covariant) = cv_2d * cv_2d * v(contravariant) |
|
|
|
|
|
! on a l'application ( x(X) , y(Y) ) avec - im/2 +1 < X < im/2 |
|
|
! = = |
|
|
! et - jm/2 < Y < jm/2 |
|
|
! = = |
|
|
|
|
|
! . x est la longitude du point en radians . |
|
|
! . y est la latitude du point en radians . |
|
|
! . . |
|
|
! . on a : cu_2d(i, j) = rad * COS(y) * dx/dX . |
|
|
! . cv( j ) = rad * dy/dY . |
|
|
! . aire_2d(i, j) = cu_2d(i, j) * cv(j) . |
|
|
! . . |
|
|
! . y, dx/dX, dy/dY calcules aux points concernes . |
|
|
! , |
|
|
! cv , bien que dependant de j uniquement, sera ici indice aussi en i |
|
|
! pour un adressage plus facile en ij . |
|
|
|
|
|
! ************** aux points u et v , ***************** |
|
|
! xprimu et xprimv sont respectivement les valeurs de dx/dX |
|
|
! yprimu et yprimv . . . . . . . . . . . dy/dY |
|
|
! rlatu et rlatv . . . . . . . . . . .la latitude |
|
|
! cvu et cv_2d . . . . . . . . . . . cv_2d |
|
|
|
|
|
! ************** aux points u, v, scalaires, et z **************** |
|
|
! cu_2d, cuv, cuscal, cuz sont respectiv. les valeurs de cu_2d |
|
|
|
|
|
! 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) |
|
|
|
|
|
USE dimens_m, ONLY : iim, jjm |
|
|
USE paramet_m, ONLY : iip1, jjp1 |
|
|
USE comconst, ONLY : g, omeg, pi, rad |
|
|
USE comdissnew, ONLY : coefdis, nitergdiv, nitergrot, niterh |
|
|
USE logic, ONLY : fxyhypb, ysinus |
|
55 |
USE comgeom, ONLY : airesurg_2d, aireu_2d, airev_2d, aire_2d, & |
USE comgeom, ONLY : airesurg_2d, aireu_2d, airev_2d, aire_2d, & |
56 |
alpha1p2_2d, alpha1p4_2d, alpha1_2d, & |
alpha1p2_2d, alpha1p4_2d, alpha1_2d, & |
57 |
alpha2p3_2d, alpha2_2d, alpha3p4_2d, alpha3_2d, alpha4_2d, apoln, & |
alpha2p3_2d, alpha2_2d, alpha3p4_2d, alpha3_2d, alpha4_2d, apoln, & |
61 |
rlatv, rlonu, rlonv, unsairez_2d, unsaire_2d, unsairz_gam_2d, & |
rlatv, rlonu, rlonv, unsairez_2d, unsaire_2d, unsairz_gam_2d, & |
62 |
unsair_gam1_2d, unsair_gam2_2d, unsapolnga1, unsapolnga2, & |
unsair_gam1_2d, unsair_gam2_2d, unsapolnga1, unsapolnga2, & |
63 |
unsapolsga1, unsapolsga2, unscu2_2d, unscv2_2d, xprimu, xprimv |
unsapolsga1, unsapolsga2, unscu2_2d, unscv2_2d, xprimu, xprimv |
64 |
|
USE comdissnew, ONLY : coefdis, nitergdiv, nitergrot, niterh |
65 |
|
USE dimens_m, ONLY : iim, jjm |
66 |
|
USE logic, ONLY : fxyhypb, ysinus |
67 |
|
use nr_util, only: pi |
68 |
|
USE paramet_m, ONLY : iip1, jjp1 |
69 |
USE serre, ONLY : alphax, alphay, clat, clon, dzoomx, dzoomy, grossismx, & |
USE serre, ONLY : alphax, alphay, clat, clon, dzoomx, dzoomy, grossismx, & |
70 |
grossismy, pxo, pyo, taux, tauy, transx, transy |
grossismy, pxo, pyo, taux, tauy, transx, transy |
71 |
|
|
93 |
real aireij2_2d(iim + 1, jjm + 1) |
real aireij2_2d(iim + 1, jjm + 1) |
94 |
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) |
95 |
real airuscv2_2d(iim + 1, jjm) |
real airuscv2_2d(iim + 1, jjm) |
96 |
real airvscu2_2d(iim + 1, jjm), aiuscv2gam_2d(iim + 1, jjm) |
real airvscu2_2d(iim + 1, jjm), aiuscv2gam_2d(iim + 1, jjm) |
97 |
real aivscu2gam_2d(iim + 1, jjm) |
real aivscu2gam_2d(iim + 1, jjm) |
98 |
|
|
99 |
!------------------------------------------------------------------ |
!------------------------------------------------------------------ |
122 |
|
|
123 |
WRITE (6, 990) |
WRITE (6, 990) |
124 |
|
|
125 |
IF ( .NOT. fxyhypb) THEN |
IF (.NOT. fxyhypb) THEN |
126 |
IF (ysinus) THEN |
IF (ysinus) THEN |
127 |
print *, ' *** Inigeom , Y = Sinus ( Latitude ) *** ' |
print *, ' Inigeom, Y = Sinus (Latitude) ' |
128 |
|
! utilisation de f(x, y) avec y = sinus de la latitude |
|
! utilisation de f(x, y ) avec y = sinus de la latitude ... |
|
|
|
|
129 |
CALL fxysinus(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, & |
CALL fxysinus(rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, & |
130 |
rlatu2, yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, & |
rlatu2, yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, & |
131 |
xprimm025, rlonp025, xprimp025) |
xprimm025, rlonp025, xprimp025) |
132 |
ELSE |
ELSE |
133 |
print *, '*** Inigeom , Y = Latitude , der. sinusoid . ***' |
print *, 'Inigeom, Y = Latitude, der. sinusoid .' |
134 |
! utilisation de f(x, y) a tangente sinusoidale , y etant la latit |
! utilisation de f(x, y) a tangente sinusoidale, y etant la latit |
135 |
|
|
136 |
pxo = clon*pi/180. |
pxo = clon*pi/180. |
137 |
pyo = 2.*clat*pi/180. |
pyo = 2.*clat*pi/180. |
138 |
|
|
139 |
! determination de transx ( pour le zoom ) par Newton-Raphson . |
! determination de transx (pour le zoom) par Newton-Raphson |
140 |
|
|
141 |
itmax = 10 |
itmax = 10 |
142 |
eps = .1E-7 |
eps = .1E-7 |
175 |
rlonp025, xprimp025) |
rlonp025, xprimp025) |
176 |
END IF |
END IF |
177 |
ELSE |
ELSE |
178 |
! .... Utilisation de fxyhyper , f(x, y) a derivee tangente hyperbol. |
! Utilisation de fxyhyper, f(x, y) à dérivée tangente hyperbolique |
179 |
print *, '*** Inigeom , Y = Latitude , der.tg. hyperbolique ***' |
print *, 'Inigeom, Y = Latitude, dérivée tangente hyperbolique' |
180 |
CALL fxyhyper(clat, grossismy, dzoomy, tauy, clon, grossismx, dzoomx, & |
CALL fxyhyper(clat, grossismy, dzoomy, tauy, clon, grossismx, dzoomx, & |
181 |
taux, rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, & |
taux, rlatu, yprimu, rlatv, yprimv, rlatu1, yprimu1, rlatu2, & |
182 |
yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, & |
yprimu2, rlonu, xprimu, rlonv, xprimv, rlonm025, xprimm025, & |
186 |
rlatu(1) = asin(1.) |
rlatu(1) = asin(1.) |
187 |
rlatu(jjp1) = -rlatu(1) |
rlatu(jjp1) = -rlatu(1) |
188 |
|
|
189 |
! .... calcul aux poles .... |
! calcul aux poles |
190 |
|
|
191 |
yprimu(1) = 0. |
yprimu(1) = 0. |
192 |
yprimu(jjp1) = 0. |
yprimu(jjp1) = 0. |
193 |
|
|
194 |
un4rad2 = 0.25*rad*rad |
un4rad2 = 0.25*rad*rad |
195 |
|
|
196 |
! 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) |
197 |
! - et de fext_2d , force de coriolis extensive . |
! - et de fext_2d, force de coriolis extensive |
198 |
|
|
199 |
! 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 |
200 |
! affectees 4 aires entourant P , calculees respectivement aux points |
! affectees 4 aires entourant P, calculees respectivement aux points |
201 |
! ( i + 1/4, j - 1/4 ) : aireij1_2d (i, j) |
! (i + 1/4, j - 1/4) : aireij1_2d (i, j) |
202 |
! ( i + 1/4, j + 1/4 ) : aireij2_2d (i, j) |
! (i + 1/4, j + 1/4) : aireij2_2d (i, j) |
203 |
! ( i - 1/4, j + 1/4 ) : aireij3_2d (i, j) |
! (i - 1/4, j + 1/4) : aireij3_2d (i, j) |
204 |
! ( i - 1/4, j - 1/4 ) : aireij4_2d (i, j) |
! (i - 1/4, j - 1/4) : aireij4_2d (i, j) |
205 |
|
|
206 |
! , |
!, |
207 |
! 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). |
208 |
! 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 |
209 |
! 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 |
210 |
! affectees au |
! affectees au |
211 |
! point (i, j) . |
! point (i, j). |
212 |
! On definit en outre les coefficients alpha comme etant egaux a |
! On definit en outre les coefficients alpha comme etant egaux a |
213 |
! (aireij / aire_2d), c.a.d par exp. |
! (aireij / aire_2d), c.a.d par exp. |
214 |
! alpha1_2d(i, j)=aireij1_2d(i, j)/aire_2d(i, j) |
! alpha1_2d(i, j)=aireij1_2d(i, j)/aire_2d(i, j) |
215 |
|
|
216 |
! 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 |
217 |
! 4 aires aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d entourant |
! 4 aires aireij1_2d, aireij2_2d, aireij3_2d, aireij4_2d entourant |
218 |
! le point U. |
! le point U. |
219 |
! Idem pour airev_2d, airez . |
! Idem pour airev_2d, airez. |
220 |
|
|
221 |
! On a , pour chaque maille : dX = dY = 1 |
! On a, pour chaque maille : dX = dY = 1 |
222 |
|
|
223 |
! . V |
! V |
224 |
|
|
225 |
! aireij4_2d . . aireij1_2d |
! aireij4_2d . . aireij1_2d |
226 |
|
|
227 |
! U . . P . U |
! U . . P . U |
228 |
|
|
229 |
! aireij3_2d . . aireij2_2d |
! aireij3_2d . . aireij2_2d |
230 |
|
|
231 |
! . V |
! V |
232 |
|
|
233 |
! Calcul des 4 aires elementaires aireij1_2d, aireij2_2d, |
! Calcul des 4 aires elementaires aireij1_2d, aireij2_2d, |
234 |
! aireij3_2d, aireij4_2d |
! aireij3_2d, aireij4_2d |
235 |
! qui entourent chaque aire_2d(i, j) , ainsi que les 4 elongations |
! qui entourent chaque aire_2d(i, j), ainsi que les 4 elongations |
236 |
! elementaires |
! elementaires |
237 |
! cuij et les 4 elongat. cvij qui sont calculees aux memes |
! cuij et les 4 elongat. cvij qui sont calculees aux memes |
238 |
! endroits que les aireij . |
! endroits que les aireij. |
239 |
|
|
240 |
! ....... do 35 : boucle sur les jjm + 1 latitudes ..... |
! do 35 : boucle sur les jjm + 1 latitudes |
241 |
|
|
242 |
DO j = 1, jjp1 |
DO j = 1, jjp1 |
243 |
|
|
334 |
|
|
335 |
END IF |
END IF |
336 |
|
|
337 |
! ........ periodicite ............ |
! periodicite |
338 |
|
|
339 |
cvij1(iip1, j) = cvij1(1, j) |
cvij1(iip1, j) = cvij1(1, j) |
340 |
cvij2(iip1, j) = cvij2(1, j) |
cvij2(iip1, j) = cvij2(1, j) |
412 |
|
|
413 |
END DO |
END DO |
414 |
|
|
415 |
! ..... Calcul des elongations cu_2d, cv_2d, cvu ......... |
! Calcul des elongations cu_2d, cv_2d, cvu |
416 |
|
|
417 |
DO j = 1, jjm |
DO j = 1, jjm |
418 |
DO i = 1, iim |
DO i = 1, iim |
460 |
cuscvugam_2d(iip1, j) = cuscvugam_2d(1, j) |
cuscvugam_2d(iip1, j) = cuscvugam_2d(1, j) |
461 |
END DO |
END DO |
462 |
|
|
463 |
! .... calcul aux poles .... |
! calcul aux poles |
464 |
|
|
465 |
DO i = 1, iip1 |
DO i = 1, iip1 |
466 |
cu_2d(i, 1) = 0. |
cu_2d(i, 1) = 0. |
490 |
aiuscv2gam_2d(iip1, j) = aiuscv2gam_2d(1, j) |
aiuscv2gam_2d(iip1, j) = aiuscv2gam_2d(1, j) |
491 |
END DO |
END DO |
492 |
|
|
493 |
! calcul des aires aux poles : |
! calcul des aires aux poles : |
494 |
|
|
495 |
apoln = sum(aire_2d(:iim, 1)) |
apoln = sum(aire_2d(:iim, 1)) |
496 |
apols = sum(aire_2d(:iim, jjp1)) |
apols = sum(aire_2d(:iim, jjp1)) |
499 |
unsapolnga2 = 1./(apoln**(-gamdi_h)) |
unsapolnga2 = 1./(apoln**(-gamdi_h)) |
500 |
unsapolsga2 = 1./(apols**(-gamdi_h)) |
unsapolsga2 = 1./(apols**(-gamdi_h)) |
501 |
|
|
502 |
! changement F. Hourdin calcul conservatif pour fext_2d |
! changement F. Hourdin calcul conservatif pour fext_2d |
503 |
! constang_2d contient le produit a * cos ( latitude ) * omega |
! constang_2d contient le produit a * cos (latitude) * omega |
504 |
|
|
505 |
DO i = 1, iim |
DO i = 1, iim |
506 |
constang_2d(i, 1) = 0. |
constang_2d(i, 1) = 0. |
514 |
constang_2d(i, jjp1) = 0. |
constang_2d(i, jjp1) = 0. |
515 |
END DO |
END DO |
516 |
|
|
517 |
! periodicite en longitude |
! periodicite en longitude |
518 |
|
|
519 |
DO j = 1, jjm |
DO j = 1, jjm |
520 |
fext_2d(iip1, j) = fext_2d(1, j) |
fext_2d(iip1, j) = fext_2d(1, j) |
525 |
|
|
526 |
! fin du changement |
! fin du changement |
527 |
|
|
528 |
print *, ' *** Coordonnees de la grille *** ' |
print *, ' Coordonnees de la grille ' |
529 |
WRITE (6, 995) |
WRITE (6, 995) |
530 |
|
|
531 |
print *, ' LONGITUDES aux pts. V ( degres ) ' |
print *, ' LONGITUDES aux pts. V (degres) ' |
532 |
WRITE (6, 995) |
WRITE (6, 995) |
533 |
DO i = 1, iip1 |
DO i = 1, iip1 |
534 |
rlonvv(i) = rlonv(i)*180./pi |
rlonvv(i) = rlonv(i)*180./pi |
536 |
WRITE (6, 400) rlonvv |
WRITE (6, 400) rlonvv |
537 |
|
|
538 |
WRITE (6, 995) |
WRITE (6, 995) |
539 |
print *, ' LATITUDES aux pts. V ( degres ) ' |
print *, ' LATITUDES aux pts. V (degres) ' |
540 |
WRITE (6, 995) |
WRITE (6, 995) |
541 |
DO i = 1, jjm |
DO i = 1, jjm |
542 |
rlatuu(i) = rlatv(i)*180./pi |
rlatuu(i) = rlatv(i)*180./pi |
547 |
rlonvv(i) = rlonu(i)*180./pi |
rlonvv(i) = rlonu(i)*180./pi |
548 |
END DO |
END DO |
549 |
WRITE (6, 995) |
WRITE (6, 995) |
550 |
print *, ' LONGITUDES aux pts. U ( degres ) ' |
print *, ' LONGITUDES aux pts. U (degres) ' |
551 |
WRITE (6, 995) |
WRITE (6, 995) |
552 |
WRITE (6, 400) rlonvv |
WRITE (6, 400) rlonvv |
553 |
WRITE (6, 995) |
WRITE (6, 995) |
554 |
|
|
555 |
print *, ' LATITUDES aux pts. U ( degres ) ' |
print *, ' LATITUDES aux pts. U (degres) ' |
556 |
WRITE (6, 995) |
WRITE (6, 995) |
557 |
DO i = 1, jjp1 |
DO i = 1, jjp1 |
558 |
rlatuu(i) = rlatu(i)*180./pi |
rlatuu(i) = rlatu(i)*180./pi |