| 1 | ! |
|---|
| 2 | ! $Id: inigeom.F 2603 2016-07-25 09:31:56Z snguyen $ |
|---|
| 3 | ! |
|---|
| 4 | c |
|---|
| 5 | c |
|---|
| 6 | SUBROUTINE inigeom |
|---|
| 7 | c |
|---|
| 8 | c Auteur : P. Le Van |
|---|
| 9 | c |
|---|
| 10 | c ............ Version du 01/04/2001 ........................ |
|---|
| 11 | c |
|---|
| 12 | c Calcul des elongations cuij1,.cuij4 , cvij1,..cvij4 aux memes en- |
|---|
| 13 | c endroits que les aires aireij1,..aireij4 . |
|---|
| 14 | |
|---|
| 15 | c Choix entre f(y) a derivee sinusoid. ou a derivee tangente hyperbol. |
|---|
| 16 | c |
|---|
| 17 | c |
|---|
| 18 | use fxhyp_m, only: fxhyp |
|---|
| 19 | use fyhyp_m, only: fyhyp |
|---|
| 20 | USE comconst_mod, ONLY: pi, g, omeg, rad |
|---|
| 21 | USE logic_mod, ONLY: fxyhypb, ysinus |
|---|
| 22 | USE serre_mod, ONLY: clon,clat,grossismx,grossismy,dzoomx,dzoomy, |
|---|
| 23 | & alphax,alphay,taux,tauy,transx,transy,pxo,pyo |
|---|
| 24 | IMPLICIT NONE |
|---|
| 25 | c |
|---|
| 26 | include "dimensions.h" |
|---|
| 27 | include "paramet.h" |
|---|
| 28 | include "comgeom2.h" |
|---|
| 29 | include "comdissnew.h" |
|---|
| 30 | |
|---|
| 31 | c----------------------------------------------------------------------- |
|---|
| 32 | c .... Variables locales .... |
|---|
| 33 | c |
|---|
| 34 | INTEGER i,j,itmax,itmay,iter |
|---|
| 35 | REAL cvu(iip1,jjp1),cuv(iip1,jjm) |
|---|
| 36 | REAL ai14,ai23,airez,rlatp,rlatm,xprm,xprp,un4rad2,yprp,yprm |
|---|
| 37 | REAL eps,x1,xo1,f,df,xdm,y1,yo1,ydm |
|---|
| 38 | REAL coslatm,coslatp,radclatm,radclatp |
|---|
| 39 | REAL cuij1(iip1,jjp1),cuij2(iip1,jjp1),cuij3(iip1,jjp1), |
|---|
| 40 | * cuij4(iip1,jjp1) |
|---|
| 41 | REAL cvij1(iip1,jjp1),cvij2(iip1,jjp1),cvij3(iip1,jjp1), |
|---|
| 42 | * cvij4(iip1,jjp1) |
|---|
| 43 | REAL rlonvv(iip1),rlatuu(jjp1) |
|---|
| 44 | REAL rlatu1(jjm),yprimu1(jjm),rlatu2(jjm),yprimu2(jjm) , |
|---|
| 45 | * yprimv(jjm),yprimu(jjp1) |
|---|
| 46 | REAL gamdi_gdiv, gamdi_grot, gamdi_h |
|---|
| 47 | |
|---|
| 48 | REAL rlonm025(iip1),xprimm025(iip1), rlonp025(iip1), |
|---|
| 49 | , xprimp025(iip1) |
|---|
| 50 | SAVE rlatu1,yprimu1,rlatu2,yprimu2,yprimv,yprimu |
|---|
| 51 | SAVE rlonm025,xprimm025,rlonp025,xprimp025 |
|---|
| 52 | |
|---|
| 53 | REAL SSUM |
|---|
| 54 | c |
|---|
| 55 | c |
|---|
| 56 | c ------------------------------------------------------------------ |
|---|
| 57 | c - - |
|---|
| 58 | c - calcul des coeff. ( cu, cv , 1./cu**2, 1./cv**2 ) - |
|---|
| 59 | c - - |
|---|
| 60 | c ------------------------------------------------------------------ |
|---|
| 61 | c |
|---|
| 62 | c les coef. ( cu, cv ) permettent de passer des vitesses naturelles |
|---|
| 63 | c aux vitesses covariantes et contravariantes , ou vice-versa ... |
|---|
| 64 | c |
|---|
| 65 | c |
|---|
| 66 | c on a : u (covariant) = cu * u (naturel) , u(contrav)= u(nat)/cu |
|---|
| 67 | c v (covariant) = cv * v (naturel) , v(contrav)= v(nat)/cv |
|---|
| 68 | c |
|---|
| 69 | c on en tire : u(covariant) = cu * cu * u(contravariant) |
|---|
| 70 | c v(covariant) = cv * cv * v(contravariant) |
|---|
| 71 | c |
|---|
| 72 | c |
|---|
| 73 | c on a l'application ( x(X) , y(Y) ) avec - im/2 +1 < X < im/2 |
|---|
| 74 | c = = |
|---|
| 75 | c et - jm/2 < Y < jm/2 |
|---|
| 76 | c = = |
|---|
| 77 | c |
|---|
| 78 | c ................................................... |
|---|
| 79 | c ................................................... |
|---|
| 80 | c . x est la longitude du point en radians . |
|---|
| 81 | c . y est la latitude du point en radians . |
|---|
| 82 | c . . |
|---|
| 83 | c . on a : cu(i,j) = rad * COS(y) * dx/dX . |
|---|
| 84 | c . cv( j ) = rad * dy/dY . |
|---|
| 85 | c . aire(i,j) = cu(i,j) * cv(j) . |
|---|
| 86 | c . . |
|---|
| 87 | c . y, dx/dX, dy/dY calcules aux points concernes . |
|---|
| 88 | c . . |
|---|
| 89 | c ................................................... |
|---|
| 90 | c ................................................... |
|---|
| 91 | c |
|---|
| 92 | c |
|---|
| 93 | c |
|---|
| 94 | c , |
|---|
| 95 | c cv , bien que dependant de j uniquement,sera ici indice aussi en i |
|---|
| 96 | c pour un adressage plus facile en ij . |
|---|
| 97 | c |
|---|
| 98 | c |
|---|
| 99 | c |
|---|
| 100 | c ************** aux points u et v , ***************** |
|---|
| 101 | c xprimu et xprimv sont respectivement les valeurs de dx/dX |
|---|
| 102 | c yprimu et yprimv . . . . . . . . . . . dy/dY |
|---|
| 103 | c rlatu et rlatv . . . . . . . . . . .la latitude |
|---|
| 104 | c cvu et cv . . . . . . . . . . . cv |
|---|
| 105 | c |
|---|
| 106 | c ************** aux points u, v, scalaires, et z **************** |
|---|
| 107 | c cu, cuv, cuscal, cuz sont respectiv. les valeurs de cu |
|---|
| 108 | c |
|---|
| 109 | c |
|---|
| 110 | c |
|---|
| 111 | c Exemple de distribution de variables sur la grille dans le |
|---|
| 112 | c domaine de travail ( X,Y ) . |
|---|
| 113 | c ................................................................ |
|---|
| 114 | c DX=DY= 1 |
|---|
| 115 | c |
|---|
| 116 | c |
|---|
| 117 | c + represente un point scalaire ( p.exp la pression ) |
|---|
| 118 | c > represente la composante zonale du vent |
|---|
| 119 | c V represente la composante meridienne du vent |
|---|
| 120 | c o represente la vorticite |
|---|
| 121 | c |
|---|
| 122 | c ---- , car aux poles , les comp.zonales covariantes sont nulles |
|---|
| 123 | c |
|---|
| 124 | c |
|---|
| 125 | c |
|---|
| 126 | c i -> |
|---|
| 127 | c |
|---|
| 128 | c 1 2 3 4 5 6 7 8 |
|---|
| 129 | c j |
|---|
| 130 | c v 1 + ---- + ---- + ---- + ---- + ---- + ---- + ---- + -- |
|---|
| 131 | c |
|---|
| 132 | c V o V o V o V o V o V o V o V o |
|---|
| 133 | c |
|---|
| 134 | c 2 + > + > + > + > + > + > + > + > |
|---|
| 135 | c |
|---|
| 136 | c V o V o V o V o V o V o V o V o |
|---|
| 137 | c |
|---|
| 138 | c 3 + > + > + > + > + > + > + > + > |
|---|
| 139 | c |
|---|
| 140 | c V o V o V o V o V o V o V o V o |
|---|
| 141 | c |
|---|
| 142 | c 4 + > + > + > + > + > + > + > + > |
|---|
| 143 | c |
|---|
| 144 | c V o V o V o V o V o V o V o V o |
|---|
| 145 | c |
|---|
| 146 | c 5 + ---- + ---- + ---- + ---- + ---- + ---- + ---- + -- |
|---|
| 147 | c |
|---|
| 148 | c |
|---|
| 149 | c Ci-dessus, on voit que le nombre de pts.en longitude est egal |
|---|
| 150 | c a IM = 8 |
|---|
| 151 | c De meme , le nombre d'intervalles entre les 2 poles est egal |
|---|
| 152 | c a JM = 4 |
|---|
| 153 | c |
|---|
| 154 | c Les points scalaires ( + ) correspondent donc a des valeurs |
|---|
| 155 | c entieres de i ( 1 a IM ) et de j ( 1 a JM +1 ) . |
|---|
| 156 | c |
|---|
| 157 | c Les vents U ( > ) correspondent a des valeurs semi- |
|---|
| 158 | c entieres de i ( 1+ 0.5 a IM+ 0.5) et entieres de j ( 1 a JM+1) |
|---|
| 159 | c |
|---|
| 160 | c Les vents V ( V ) correspondent a des valeurs entieres |
|---|
| 161 | c de i ( 1 a IM ) et semi-entieres de j ( 1 +0.5 a JM +0.5) |
|---|
| 162 | c |
|---|
| 163 | c |
|---|
| 164 | c |
|---|
| 165 | WRITE(6,3) |
|---|
| 166 | 3 FORMAT( // 10x,' .... INIGEOM date du 01/06/98 ..... ', |
|---|
| 167 | * //5x,' Calcul des elongations cu et cv comme sommes des 4 ' / |
|---|
| 168 | * 5x,' elong. cuij1, .. 4 , cvij1,.. 4 qui les entourent , aux |
|---|
| 169 | * '/ 5x,' memes endroits que les aires aireij1,...j4 . ' / ) |
|---|
| 170 | c |
|---|
| 171 | c |
|---|
| 172 | IF( nitergdiv.NE.2 ) THEN |
|---|
| 173 | gamdi_gdiv = coefdis/ ( REAL(nitergdiv) -2. ) |
|---|
| 174 | ELSE |
|---|
| 175 | gamdi_gdiv = 0. |
|---|
| 176 | ENDIF |
|---|
| 177 | IF( nitergrot.NE.2 ) THEN |
|---|
| 178 | gamdi_grot = coefdis/ ( REAL(nitergrot) -2. ) |
|---|
| 179 | ELSE |
|---|
| 180 | gamdi_grot = 0. |
|---|
| 181 | ENDIF |
|---|
| 182 | IF( niterh.NE.2 ) THEN |
|---|
| 183 | gamdi_h = coefdis/ ( REAL(niterh) -2. ) |
|---|
| 184 | ELSE |
|---|
| 185 | gamdi_h = 0. |
|---|
| 186 | ENDIF |
|---|
| 187 | |
|---|
| 188 | WRITE(6,*) ' gamdi_gd ',gamdi_gdiv,gamdi_grot,gamdi_h,coefdis, |
|---|
| 189 | * nitergdiv,nitergrot,niterh |
|---|
| 190 | c |
|---|
| 191 | pi = 2.* ASIN(1.) |
|---|
| 192 | c |
|---|
| 193 | WRITE(6,990) |
|---|
| 194 | |
|---|
| 195 | c ---------------------------------------------------------------- |
|---|
| 196 | c |
|---|
| 197 | IF( .NOT.fxyhypb ) THEN |
|---|
| 198 | c |
|---|
| 199 | c |
|---|
| 200 | IF( ysinus ) THEN |
|---|
| 201 | c |
|---|
| 202 | WRITE(6,*) ' *** Inigeom , Y = Sinus ( Latitude ) *** ' |
|---|
| 203 | c |
|---|
| 204 | c .... utilisation de f(x,y ) avec y = sinus de la latitude ..... |
|---|
| 205 | |
|---|
| 206 | CALL fxysinus (rlatu,yprimu,rlatv,yprimv,rlatu1,yprimu1, |
|---|
| 207 | , rlatu2,yprimu2, |
|---|
| 208 | , rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025) |
|---|
| 209 | |
|---|
| 210 | ELSE |
|---|
| 211 | c |
|---|
| 212 | WRITE(6,*) '*** Inigeom , Y = Latitude , der. sinusoid . ***' |
|---|
| 213 | |
|---|
| 214 | c .... utilisation de f(x,y) a tangente sinusoidale , y etant la latit. ... |
|---|
| 215 | c |
|---|
| 216 | |
|---|
| 217 | pxo = clon *pi /180. |
|---|
| 218 | pyo = 2.* clat* pi /180. |
|---|
| 219 | c |
|---|
| 220 | c .... determination de transx ( pour le zoom ) par Newton-Raphson ... |
|---|
| 221 | c |
|---|
| 222 | itmax = 10 |
|---|
| 223 | eps = .1e-7 |
|---|
| 224 | c |
|---|
| 225 | xo1 = 0. |
|---|
| 226 | DO 10 iter = 1, itmax |
|---|
| 227 | x1 = xo1 |
|---|
| 228 | f = x1+ alphax *SIN(x1-pxo) |
|---|
| 229 | df = 1.+ alphax *COS(x1-pxo) |
|---|
| 230 | x1 = x1 - f/df |
|---|
| 231 | xdm = ABS( x1- xo1 ) |
|---|
| 232 | IF( xdm.LE.eps )GO TO 11 |
|---|
| 233 | xo1 = x1 |
|---|
| 234 | 10 CONTINUE |
|---|
| 235 | 11 CONTINUE |
|---|
| 236 | c |
|---|
| 237 | transx = xo1 |
|---|
| 238 | |
|---|
| 239 | itmay = 10 |
|---|
| 240 | eps = .1e-7 |
|---|
| 241 | C |
|---|
| 242 | yo1 = 0. |
|---|
| 243 | DO 15 iter = 1,itmay |
|---|
| 244 | y1 = yo1 |
|---|
| 245 | f = y1 + alphay* SIN(y1-pyo) |
|---|
| 246 | df = 1. + alphay* COS(y1-pyo) |
|---|
| 247 | y1 = y1 -f/df |
|---|
| 248 | ydm = ABS(y1-yo1) |
|---|
| 249 | IF(ydm.LE.eps) GO TO 17 |
|---|
| 250 | yo1 = y1 |
|---|
| 251 | 15 CONTINUE |
|---|
| 252 | c |
|---|
| 253 | 17 CONTINUE |
|---|
| 254 | transy = yo1 |
|---|
| 255 | |
|---|
| 256 | CALL fxy ( rlatu,yprimu,rlatv,yprimv,rlatu1,yprimu1, |
|---|
| 257 | , rlatu2,yprimu2, |
|---|
| 258 | , rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025) |
|---|
| 259 | |
|---|
| 260 | ENDIF |
|---|
| 261 | c |
|---|
| 262 | ELSE |
|---|
| 263 | c |
|---|
| 264 | c .... Utilisation de fxyhyper , f(x,y) a derivee tangente hyperbol. |
|---|
| 265 | c ..................................................................... |
|---|
| 266 | |
|---|
| 267 | WRITE(6,*)'*** Inigeom , Y = Latitude , der.tg. hyperbolique ***' |
|---|
| 268 | |
|---|
| 269 | CALL fyhyp(rlatu, yprimu, rlatv, rlatu2, yprimu2, rlatu1, yprimu1) |
|---|
| 270 | CALL fxhyp(xprimm025, rlonv, xprimv, rlonu, xprimu, xprimp025) |
|---|
| 271 | |
|---|
| 272 | ENDIF |
|---|
| 273 | c |
|---|
| 274 | c ------------------------------------------------------------------- |
|---|
| 275 | |
|---|
| 276 | c |
|---|
| 277 | rlatu(1) = ASIN(1.) |
|---|
| 278 | rlatu(jjp1) = - rlatu(1) |
|---|
| 279 | c |
|---|
| 280 | c |
|---|
| 281 | c .... calcul aux poles .... |
|---|
| 282 | c |
|---|
| 283 | yprimu(1) = 0. |
|---|
| 284 | yprimu(jjp1) = 0. |
|---|
| 285 | c |
|---|
| 286 | c |
|---|
| 287 | un4rad2 = 0.25 * rad * rad |
|---|
| 288 | c |
|---|
| 289 | c -------------------------------------------------------------------- |
|---|
| 290 | c -------------------------------------------------------------------- |
|---|
| 291 | c - - |
|---|
| 292 | c - calcul des aires ( aire,aireu,airev, 1./aire, 1./airez ) - |
|---|
| 293 | c - et de fext , force de coriolis extensive . - |
|---|
| 294 | c - - |
|---|
| 295 | c -------------------------------------------------------------------- |
|---|
| 296 | c -------------------------------------------------------------------- |
|---|
| 297 | c |
|---|
| 298 | c |
|---|
| 299 | c |
|---|
| 300 | c A 1 point scalaire P (i,j) de la grille, reguliere en (X,Y) , sont |
|---|
| 301 | c affectees 4 aires entourant P , calculees respectivement aux points |
|---|
| 302 | c ( i + 1/4, j - 1/4 ) : aireij1 (i,j) |
|---|
| 303 | c ( i + 1/4, j + 1/4 ) : aireij2 (i,j) |
|---|
| 304 | c ( i - 1/4, j + 1/4 ) : aireij3 (i,j) |
|---|
| 305 | c ( i - 1/4, j - 1/4 ) : aireij4 (i,j) |
|---|
| 306 | c |
|---|
| 307 | c , |
|---|
| 308 | c Les cotes de chacun de ces 4 carres etant egaux a 1/2 suivant (X,Y). |
|---|
| 309 | c Chaque aire centree en 1 point scalaire P(i,j) est egale a la somme |
|---|
| 310 | c des 4 aires aireij1,aireij2,aireij3,aireij4 qui sont affectees au |
|---|
| 311 | c point (i,j) . |
|---|
| 312 | c On definit en outre les coefficients alpha comme etant egaux a |
|---|
| 313 | c (aireij / aire), c.a.d par exp. alpha1(i,j)=aireij1(i,j)/aire(i,j) |
|---|
| 314 | c |
|---|
| 315 | c De meme, toute aire centree en 1 point U est egale a la somme des |
|---|
| 316 | c 4 aires aireij1,aireij2,aireij3,aireij4 entourant le point U . |
|---|
| 317 | c Idem pour airev, airez . |
|---|
| 318 | c |
|---|
| 319 | c On a ,pour chaque maille : dX = dY = 1 |
|---|
| 320 | c |
|---|
| 321 | c |
|---|
| 322 | c . V |
|---|
| 323 | c |
|---|
| 324 | c aireij4 . . aireij1 |
|---|
| 325 | c |
|---|
| 326 | c U . . P . U |
|---|
| 327 | c |
|---|
| 328 | c aireij3 . . aireij2 |
|---|
| 329 | c |
|---|
| 330 | c . V |
|---|
| 331 | c |
|---|
| 332 | c |
|---|
| 333 | c |
|---|
| 334 | c |
|---|
| 335 | c |
|---|
| 336 | c .................................................................... |
|---|
| 337 | c |
|---|
| 338 | c Calcul des 4 aires elementaires aireij1,aireij2,aireij3,aireij4 |
|---|
| 339 | c qui entourent chaque aire(i,j) , ainsi que les 4 elongations elemen |
|---|
| 340 | c taires cuij et les 4 elongat. cvij qui sont calculees aux memes |
|---|
| 341 | c endroits que les aireij . |
|---|
| 342 | c |
|---|
| 343 | c .................................................................... |
|---|
| 344 | c |
|---|
| 345 | c ....... do 35 : boucle sur les jjm + 1 latitudes ..... |
|---|
| 346 | c |
|---|
| 347 | c |
|---|
| 348 | DO 35 j = 1, jjp1 |
|---|
| 349 | c |
|---|
| 350 | IF ( j. eq. 1 ) THEN |
|---|
| 351 | c |
|---|
| 352 | yprm = yprimu1(j) |
|---|
| 353 | rlatm = rlatu1(j) |
|---|
| 354 | c |
|---|
| 355 | coslatm = COS( rlatm ) |
|---|
| 356 | radclatm = 0.5* rad * coslatm |
|---|
| 357 | c |
|---|
| 358 | DO 30 i = 1, iim |
|---|
| 359 | xprp = xprimp025( i ) |
|---|
| 360 | xprm = xprimm025( i ) |
|---|
| 361 | aireij2( i,1 ) = un4rad2 * coslatm * xprp * yprm |
|---|
| 362 | aireij3( i,1 ) = un4rad2 * coslatm * xprm * yprm |
|---|
| 363 | cuij2 ( i,1 ) = radclatm * xprp |
|---|
| 364 | cuij3 ( i,1 ) = radclatm * xprm |
|---|
| 365 | cvij2 ( i,1 ) = 0.5* rad * yprm |
|---|
| 366 | cvij3 ( i,1 ) = cvij2(i,1) |
|---|
| 367 | 30 CONTINUE |
|---|
| 368 | c |
|---|
| 369 | DO i = 1, iim |
|---|
| 370 | aireij1( i,1 ) = 0. |
|---|
| 371 | aireij4( i,1 ) = 0. |
|---|
| 372 | cuij1 ( i,1 ) = 0. |
|---|
| 373 | cuij4 ( i,1 ) = 0. |
|---|
| 374 | cvij1 ( i,1 ) = 0. |
|---|
| 375 | cvij4 ( i,1 ) = 0. |
|---|
| 376 | ENDDO |
|---|
| 377 | c |
|---|
| 378 | END IF |
|---|
| 379 | c |
|---|
| 380 | IF ( j. eq. jjp1 ) THEN |
|---|
| 381 | yprp = yprimu2(j-1) |
|---|
| 382 | rlatp = rlatu2 (j-1) |
|---|
| 383 | ccc yprp = fyprim( REAL(j) - 0.25 ) |
|---|
| 384 | ccc rlatp = fy ( REAL(j) - 0.25 ) |
|---|
| 385 | c |
|---|
| 386 | coslatp = COS( rlatp ) |
|---|
| 387 | radclatp = 0.5* rad * coslatp |
|---|
| 388 | c |
|---|
| 389 | DO 31 i = 1,iim |
|---|
| 390 | xprp = xprimp025( i ) |
|---|
| 391 | xprm = xprimm025( i ) |
|---|
| 392 | aireij1( i,jjp1 ) = un4rad2 * coslatp * xprp * yprp |
|---|
| 393 | aireij4( i,jjp1 ) = un4rad2 * coslatp * xprm * yprp |
|---|
| 394 | cuij1(i,jjp1) = radclatp * xprp |
|---|
| 395 | cuij4(i,jjp1) = radclatp * xprm |
|---|
| 396 | cvij1(i,jjp1) = 0.5 * rad* yprp |
|---|
| 397 | cvij4(i,jjp1) = cvij1(i,jjp1) |
|---|
| 398 | 31 CONTINUE |
|---|
| 399 | c |
|---|
| 400 | DO i = 1, iim |
|---|
| 401 | aireij2( i,jjp1 ) = 0. |
|---|
| 402 | aireij3( i,jjp1 ) = 0. |
|---|
| 403 | cvij2 ( i,jjp1 ) = 0. |
|---|
| 404 | cvij3 ( i,jjp1 ) = 0. |
|---|
| 405 | cuij2 ( i,jjp1 ) = 0. |
|---|
| 406 | cuij3 ( i,jjp1 ) = 0. |
|---|
| 407 | ENDDO |
|---|
| 408 | c |
|---|
| 409 | END IF |
|---|
| 410 | c |
|---|
| 411 | |
|---|
| 412 | IF ( j .gt. 1 .AND. j .lt. jjp1 ) THEN |
|---|
| 413 | c |
|---|
| 414 | rlatp = rlatu2 ( j-1 ) |
|---|
| 415 | yprp = yprimu2( j-1 ) |
|---|
| 416 | rlatm = rlatu1 ( j ) |
|---|
| 417 | yprm = yprimu1( j ) |
|---|
| 418 | cc rlatp = fy ( REAL(j) - 0.25 ) |
|---|
| 419 | cc yprp = fyprim( REAL(j) - 0.25 ) |
|---|
| 420 | cc rlatm = fy ( REAL(j) + 0.25 ) |
|---|
| 421 | cc yprm = fyprim( REAL(j) + 0.25 ) |
|---|
| 422 | |
|---|
| 423 | coslatm = COS( rlatm ) |
|---|
| 424 | coslatp = COS( rlatp ) |
|---|
| 425 | radclatp = 0.5* rad * coslatp |
|---|
| 426 | radclatm = 0.5* rad * coslatm |
|---|
| 427 | c |
|---|
| 428 | ai14 = un4rad2 * coslatp * yprp |
|---|
| 429 | ai23 = un4rad2 * coslatm * yprm |
|---|
| 430 | DO 32 i = 1,iim |
|---|
| 431 | xprp = xprimp025( i ) |
|---|
| 432 | xprm = xprimm025( i ) |
|---|
| 433 | |
|---|
| 434 | aireij1 ( i,j ) = ai14 * xprp |
|---|
| 435 | aireij2 ( i,j ) = ai23 * xprp |
|---|
| 436 | aireij3 ( i,j ) = ai23 * xprm |
|---|
| 437 | aireij4 ( i,j ) = ai14 * xprm |
|---|
| 438 | cuij1 ( i,j ) = radclatp * xprp |
|---|
| 439 | cuij2 ( i,j ) = radclatm * xprp |
|---|
| 440 | cuij3 ( i,j ) = radclatm * xprm |
|---|
| 441 | cuij4 ( i,j ) = radclatp * xprm |
|---|
| 442 | cvij1 ( i,j ) = 0.5* rad * yprp |
|---|
| 443 | cvij2 ( i,j ) = 0.5* rad * yprm |
|---|
| 444 | cvij3 ( i,j ) = cvij2(i,j) |
|---|
| 445 | cvij4 ( i,j ) = cvij1(i,j) |
|---|
| 446 | 32 CONTINUE |
|---|
| 447 | c |
|---|
| 448 | END IF |
|---|
| 449 | c |
|---|
| 450 | c ........ periodicite ............ |
|---|
| 451 | c |
|---|
| 452 | cvij1 (iip1,j) = cvij1 (1,j) |
|---|
| 453 | cvij2 (iip1,j) = cvij2 (1,j) |
|---|
| 454 | cvij3 (iip1,j) = cvij3 (1,j) |
|---|
| 455 | cvij4 (iip1,j) = cvij4 (1,j) |
|---|
| 456 | cuij1 (iip1,j) = cuij1 (1,j) |
|---|
| 457 | cuij2 (iip1,j) = cuij2 (1,j) |
|---|
| 458 | cuij3 (iip1,j) = cuij3 (1,j) |
|---|
| 459 | cuij4 (iip1,j) = cuij4 (1,j) |
|---|
| 460 | aireij1 (iip1,j) = aireij1 (1,j ) |
|---|
| 461 | aireij2 (iip1,j) = aireij2 (1,j ) |
|---|
| 462 | aireij3 (iip1,j) = aireij3 (1,j ) |
|---|
| 463 | aireij4 (iip1,j) = aireij4 (1,j ) |
|---|
| 464 | |
|---|
| 465 | 35 CONTINUE |
|---|
| 466 | c |
|---|
| 467 | c .............................................................. |
|---|
| 468 | c |
|---|
| 469 | DO 37 j = 1, jjp1 |
|---|
| 470 | DO 36 i = 1, iim |
|---|
| 471 | aire ( i,j ) = aireij1(i,j) + aireij2(i,j) + aireij3(i,j) + |
|---|
| 472 | * aireij4(i,j) |
|---|
| 473 | alpha1 ( i,j ) = aireij1(i,j) / aire(i,j) |
|---|
| 474 | alpha2 ( i,j ) = aireij2(i,j) / aire(i,j) |
|---|
| 475 | alpha3 ( i,j ) = aireij3(i,j) / aire(i,j) |
|---|
| 476 | alpha4 ( i,j ) = aireij4(i,j) / aire(i,j) |
|---|
| 477 | alpha1p2( i,j ) = alpha1 (i,j) + alpha2 (i,j) |
|---|
| 478 | alpha1p4( i,j ) = alpha1 (i,j) + alpha4 (i,j) |
|---|
| 479 | alpha2p3( i,j ) = alpha2 (i,j) + alpha3 (i,j) |
|---|
| 480 | alpha3p4( i,j ) = alpha3 (i,j) + alpha4 (i,j) |
|---|
| 481 | 36 CONTINUE |
|---|
| 482 | c |
|---|
| 483 | c |
|---|
| 484 | aire (iip1,j) = aire (1,j) |
|---|
| 485 | alpha1 (iip1,j) = alpha1 (1,j) |
|---|
| 486 | alpha2 (iip1,j) = alpha2 (1,j) |
|---|
| 487 | alpha3 (iip1,j) = alpha3 (1,j) |
|---|
| 488 | alpha4 (iip1,j) = alpha4 (1,j) |
|---|
| 489 | alpha1p2(iip1,j) = alpha1p2(1,j) |
|---|
| 490 | alpha1p4(iip1,j) = alpha1p4(1,j) |
|---|
| 491 | alpha2p3(iip1,j) = alpha2p3(1,j) |
|---|
| 492 | alpha3p4(iip1,j) = alpha3p4(1,j) |
|---|
| 493 | 37 CONTINUE |
|---|
| 494 | c |
|---|
| 495 | |
|---|
| 496 | DO 42 j = 1,jjp1 |
|---|
| 497 | DO 41 i = 1,iim |
|---|
| 498 | aireu (i,j)= aireij1(i,j) + aireij2(i,j) + aireij4(i+1,j) + |
|---|
| 499 | * aireij3(i+1,j) |
|---|
| 500 | unsaire ( i,j)= 1./ aire(i,j) |
|---|
| 501 | unsair_gam1( i,j)= unsaire(i,j)** ( - gamdi_gdiv ) |
|---|
| 502 | unsair_gam2( i,j)= unsaire(i,j)** ( - gamdi_h ) |
|---|
| 503 | airesurg ( i,j)= aire(i,j)/ g |
|---|
| 504 | 41 CONTINUE |
|---|
| 505 | aireu (iip1,j) = aireu (1,j) |
|---|
| 506 | unsaire (iip1,j) = unsaire(1,j) |
|---|
| 507 | unsair_gam1(iip1,j) = unsair_gam1(1,j) |
|---|
| 508 | unsair_gam2(iip1,j) = unsair_gam2(1,j) |
|---|
| 509 | airesurg (iip1,j) = airesurg(1,j) |
|---|
| 510 | 42 CONTINUE |
|---|
| 511 | c |
|---|
| 512 | c |
|---|
| 513 | DO 48 j = 1,jjm |
|---|
| 514 | c |
|---|
| 515 | DO i=1,iim |
|---|
| 516 | airev (i,j) = aireij2(i,j)+ aireij3(i,j)+ aireij1(i,j+1) + |
|---|
| 517 | * aireij4(i,j+1) |
|---|
| 518 | ENDDO |
|---|
| 519 | DO i=1,iim |
|---|
| 520 | airez = aireij2(i,j)+aireij1(i,j+1)+aireij3(i+1,j) + |
|---|
| 521 | * aireij4(i+1,j+1) |
|---|
| 522 | unsairez(i,j) = 1./ airez |
|---|
| 523 | unsairz_gam(i,j)= unsairez(i,j)** ( - gamdi_grot ) |
|---|
| 524 | fext (i,j) = airez * SIN(rlatv(j))* 2.* omeg |
|---|
| 525 | ENDDO |
|---|
| 526 | airev (iip1,j) = airev(1,j) |
|---|
| 527 | unsairez (iip1,j) = unsairez(1,j) |
|---|
| 528 | fext (iip1,j) = fext(1,j) |
|---|
| 529 | unsairz_gam(iip1,j) = unsairz_gam(1,j) |
|---|
| 530 | c |
|---|
| 531 | 48 CONTINUE |
|---|
| 532 | c |
|---|
| 533 | c |
|---|
| 534 | c ..... Calcul des elongations cu,cv, cvu ......... |
|---|
| 535 | c |
|---|
| 536 | DO j = 1, jjm |
|---|
| 537 | DO i = 1, iim |
|---|
| 538 | cv(i,j) = 0.5 *( cvij2(i,j)+cvij3(i,j)+cvij1(i,j+1)+cvij4(i,j+1)) |
|---|
| 539 | cvu(i,j)= 0.5 *( cvij1(i,j)+cvij4(i,j)+cvij2(i,j) +cvij3(i,j) ) |
|---|
| 540 | cuv(i,j)= 0.5 *( cuij2(i,j)+cuij3(i,j)+cuij1(i,j+1)+cuij4(i,j+1)) |
|---|
| 541 | unscv2(i,j) = 1./ ( cv(i,j)*cv(i,j) ) |
|---|
| 542 | ENDDO |
|---|
| 543 | DO i = 1, iim |
|---|
| 544 | cuvsurcv (i,j) = airev(i,j) * unscv2(i,j) |
|---|
| 545 | cvsurcuv (i,j) = 1./cuvsurcv(i,j) |
|---|
| 546 | cuvscvgam1(i,j) = cuvsurcv (i,j) ** ( - gamdi_gdiv ) |
|---|
| 547 | cuvscvgam2(i,j) = cuvsurcv (i,j) ** ( - gamdi_h ) |
|---|
| 548 | cvscuvgam(i,j) = cvsurcuv (i,j) ** ( - gamdi_grot ) |
|---|
| 549 | ENDDO |
|---|
| 550 | cv (iip1,j) = cv (1,j) |
|---|
| 551 | cvu (iip1,j) = cvu (1,j) |
|---|
| 552 | unscv2 (iip1,j) = unscv2 (1,j) |
|---|
| 553 | cuv (iip1,j) = cuv (1,j) |
|---|
| 554 | cuvsurcv (iip1,j) = cuvsurcv (1,j) |
|---|
| 555 | cvsurcuv (iip1,j) = cvsurcuv (1,j) |
|---|
| 556 | cuvscvgam1(iip1,j) = cuvscvgam1(1,j) |
|---|
| 557 | cuvscvgam2(iip1,j) = cuvscvgam2(1,j) |
|---|
| 558 | cvscuvgam(iip1,j) = cvscuvgam(1,j) |
|---|
| 559 | ENDDO |
|---|
| 560 | |
|---|
| 561 | DO j = 2, jjm |
|---|
| 562 | DO i = 1, iim |
|---|
| 563 | cu(i,j) = 0.5*(cuij1(i,j)+cuij4(i+1,j)+cuij2(i,j)+cuij3(i+1,j)) |
|---|
| 564 | unscu2 (i,j) = 1./ ( cu(i,j) * cu(i,j) ) |
|---|
| 565 | cvusurcu (i,j) = aireu(i,j) * unscu2(i,j) |
|---|
| 566 | cusurcvu (i,j) = 1./ cvusurcu(i,j) |
|---|
| 567 | cvuscugam1 (i,j) = cvusurcu(i,j) ** ( - gamdi_gdiv ) |
|---|
| 568 | cvuscugam2 (i,j) = cvusurcu(i,j) ** ( - gamdi_h ) |
|---|
| 569 | cuscvugam (i,j) = cusurcvu(i,j) ** ( - gamdi_grot ) |
|---|
| 570 | ENDDO |
|---|
| 571 | cu (iip1,j) = cu(1,j) |
|---|
| 572 | unscu2 (iip1,j) = unscu2(1,j) |
|---|
| 573 | cvusurcu (iip1,j) = cvusurcu(1,j) |
|---|
| 574 | cusurcvu (iip1,j) = cusurcvu(1,j) |
|---|
| 575 | cvuscugam1(iip1,j) = cvuscugam1(1,j) |
|---|
| 576 | cvuscugam2(iip1,j) = cvuscugam2(1,j) |
|---|
| 577 | cuscvugam (iip1,j) = cuscvugam(1,j) |
|---|
| 578 | ENDDO |
|---|
| 579 | |
|---|
| 580 | c |
|---|
| 581 | c .... calcul aux poles .... |
|---|
| 582 | c |
|---|
| 583 | DO i = 1, iip1 |
|---|
| 584 | cu ( i, 1 ) = 0. |
|---|
| 585 | unscu2( i, 1 ) = 0. |
|---|
| 586 | cvu ( i, 1 ) = 0. |
|---|
| 587 | c |
|---|
| 588 | cu (i, jjp1) = 0. |
|---|
| 589 | unscu2(i, jjp1) = 0. |
|---|
| 590 | cvu (i, jjp1) = 0. |
|---|
| 591 | ENDDO |
|---|
| 592 | c |
|---|
| 593 | c .............................................................. |
|---|
| 594 | c |
|---|
| 595 | DO j = 1, jjm |
|---|
| 596 | DO i= 1, iim |
|---|
| 597 | airvscu2 (i,j) = airev(i,j)/ ( cuv(i,j) * cuv(i,j) ) |
|---|
| 598 | aivscu2gam(i,j) = airvscu2(i,j)** ( - gamdi_grot ) |
|---|
| 599 | ENDDO |
|---|
| 600 | airvscu2 (iip1,j) = airvscu2(1,j) |
|---|
| 601 | aivscu2gam(iip1,j) = aivscu2gam(1,j) |
|---|
| 602 | ENDDO |
|---|
| 603 | |
|---|
| 604 | DO j=2,jjm |
|---|
| 605 | DO i=1,iim |
|---|
| 606 | airuscv2 (i,j) = aireu(i,j)/ ( cvu(i,j) * cvu(i,j) ) |
|---|
| 607 | aiuscv2gam (i,j) = airuscv2(i,j)** ( - gamdi_grot ) |
|---|
| 608 | ENDDO |
|---|
| 609 | airuscv2 (iip1,j) = airuscv2 (1,j) |
|---|
| 610 | aiuscv2gam(iip1,j) = aiuscv2gam(1,j) |
|---|
| 611 | ENDDO |
|---|
| 612 | |
|---|
| 613 | c |
|---|
| 614 | c calcul des aires aux poles : |
|---|
| 615 | c ----------------------------- |
|---|
| 616 | c |
|---|
| 617 | apoln = SSUM(iim,aire(1,1),1) |
|---|
| 618 | apols = SSUM(iim,aire(1,jjp1),1) |
|---|
| 619 | unsapolnga1 = 1./ ( apoln ** ( - gamdi_gdiv ) ) |
|---|
| 620 | unsapolsga1 = 1./ ( apols ** ( - gamdi_gdiv ) ) |
|---|
| 621 | unsapolnga2 = 1./ ( apoln ** ( - gamdi_h ) ) |
|---|
| 622 | unsapolsga2 = 1./ ( apols ** ( - gamdi_h ) ) |
|---|
| 623 | c |
|---|
| 624 | c----------------------------------------------------------------------- |
|---|
| 625 | c gtitre='Coriolis version ancienne' |
|---|
| 626 | c gfichier='fext1' |
|---|
| 627 | c CALL writestd(fext,iip1*jjm) |
|---|
| 628 | c |
|---|
| 629 | c changement F. Hourdin calcul conservatif pour fext |
|---|
| 630 | c constang contient le produit a * cos ( latitude ) * omega |
|---|
| 631 | c |
|---|
| 632 | DO i=1,iim |
|---|
| 633 | constang(i,1) = 0. |
|---|
| 634 | ENDDO |
|---|
| 635 | DO j=1,jjm-1 |
|---|
| 636 | DO i=1,iim |
|---|
| 637 | constang(i,j+1) = rad*omeg*cu(i,j+1)*COS(rlatu(j+1)) |
|---|
| 638 | ENDDO |
|---|
| 639 | ENDDO |
|---|
| 640 | DO i=1,iim |
|---|
| 641 | constang(i,jjp1) = 0. |
|---|
| 642 | ENDDO |
|---|
| 643 | c |
|---|
| 644 | c periodicite en longitude |
|---|
| 645 | c |
|---|
| 646 | DO j=1,jjm |
|---|
| 647 | fext(iip1,j) = fext(1,j) |
|---|
| 648 | ENDDO |
|---|
| 649 | DO j=1,jjp1 |
|---|
| 650 | constang(iip1,j) = constang(1,j) |
|---|
| 651 | ENDDO |
|---|
| 652 | |
|---|
| 653 | c fin du changement |
|---|
| 654 | |
|---|
| 655 | c |
|---|
| 656 | c----------------------------------------------------------------------- |
|---|
| 657 | c |
|---|
| 658 | WRITE(6,*) ' *** Coordonnees de la grille *** ' |
|---|
| 659 | WRITE(6,995) |
|---|
| 660 | c |
|---|
| 661 | WRITE(6,*) ' LONGITUDES aux pts. V ( degres ) ' |
|---|
| 662 | WRITE(6,995) |
|---|
| 663 | DO i=1,iip1 |
|---|
| 664 | rlonvv(i) = rlonv(i)*180./pi |
|---|
| 665 | ENDDO |
|---|
| 666 | WRITE(6,400) rlonvv |
|---|
| 667 | c |
|---|
| 668 | WRITE(6,995) |
|---|
| 669 | WRITE(6,*) ' LATITUDES aux pts. V ( degres ) ' |
|---|
| 670 | WRITE(6,995) |
|---|
| 671 | DO i=1,jjm |
|---|
| 672 | rlatuu(i)=rlatv(i)*180./pi |
|---|
| 673 | ENDDO |
|---|
| 674 | WRITE(6,400) (rlatuu(i),i=1,jjm) |
|---|
| 675 | c |
|---|
| 676 | DO i=1,iip1 |
|---|
| 677 | rlonvv(i)=rlonu(i)*180./pi |
|---|
| 678 | ENDDO |
|---|
| 679 | WRITE(6,995) |
|---|
| 680 | WRITE(6,*) ' LONGITUDES aux pts. U ( degres ) ' |
|---|
| 681 | WRITE(6,995) |
|---|
| 682 | WRITE(6,400) rlonvv |
|---|
| 683 | WRITE(6,995) |
|---|
| 684 | |
|---|
| 685 | WRITE(6,*) ' LATITUDES aux pts. U ( degres ) ' |
|---|
| 686 | WRITE(6,995) |
|---|
| 687 | DO i=1,jjp1 |
|---|
| 688 | rlatuu(i)=rlatu(i)*180./pi |
|---|
| 689 | ENDDO |
|---|
| 690 | WRITE(6,400) (rlatuu(i),i=1,jjp1) |
|---|
| 691 | WRITE(6,995) |
|---|
| 692 | c |
|---|
| 693 | 444 format(f10.3,f6.0) |
|---|
| 694 | 400 FORMAT(1x,8f8.2) |
|---|
| 695 | 990 FORMAT(//) |
|---|
| 696 | 995 FORMAT(/) |
|---|
| 697 | c |
|---|
| 698 | RETURN |
|---|
| 699 | END |
|---|
