| 1 | SUBROUTINE fxyhyp( clat, rlatu,yprimu , rlatv,yprimv , |
|---|
| 2 | * rlatu1 ,yprimu1, rlatu2,yprimu2 ) |
|---|
| 3 | c |
|---|
| 4 | IMPLICIT NONE |
|---|
| 5 | c ----------------------- |
|---|
| 6 | c Auteur : P. Le Van |
|---|
| 7 | c ----------------------- |
|---|
| 8 | c |
|---|
| 9 | c .. Objet : Calcul d'une fonction f(y) ayant pour derivee une tangente |
|---|
| 10 | c hyperbolique , donc des ordonnees y de la grille horizontale , |
|---|
| 11 | c ainsi que ses derivees yprim , utlisees dans Inigeom . |
|---|
| 12 | c |
|---|
| 13 | c |
|---|
| 14 | c .... rlatu = fxyhyp ( j ) avec 1 <= j <= jjm + 1 |
|---|
| 15 | c .... rlatv = fxyhyp ( j + 0.5 ) avec 1 <= j <= jjm |
|---|
| 16 | c .... rlatu1 = fxyhyp ( j + 0.25 ) avec 1 <= j <= jjm |
|---|
| 17 | c .... rlatu2 = fxyhyp ( j + 0.75 ) avec 1 <= j <= jjm |
|---|
| 18 | c |
|---|
| 19 | |
|---|
| 20 | INTEGER nmax |
|---|
| 21 | PARAMETER ( nmax = 100 000 ) |
|---|
| 22 | c |
|---|
| 23 | #include "dimensions.h" |
|---|
| 24 | #include "paramet.h" |
|---|
| 25 | #include "comconst.h" |
|---|
| 26 | c |
|---|
| 27 | c ................. Arguments ................................ |
|---|
| 28 | c |
|---|
| 29 | c clat est (en degres) la latitude du centre du zoom : arg. d'entree |
|---|
| 30 | c Les autres arguments sont des argum. de sortie |
|---|
| 31 | c |
|---|
| 32 | REAL clat |
|---|
| 33 | REAL rlatu(jjp1),rlatv(jjm),yprimu(jjp1),yprimv(jjm) |
|---|
| 34 | REAL rlatu1(jjm),yprimu1(jjm),rlatu2(jjm),yprimu2(jjm) |
|---|
| 35 | c |
|---|
| 36 | c ................................................................. |
|---|
| 37 | c |
|---|
| 38 | c ...... . Variables locales ......... |
|---|
| 39 | c |
|---|
| 40 | REAL yprim( 0:nmax ),yi( 0:nmax ),yf( 0:nmax ) |
|---|
| 41 | REAL yyy(0:nmax) |
|---|
| 42 | REAL delta,den,eps,phij,pijjm,pis2,pis2pyo,pism,psi,psi0 |
|---|
| 43 | REAL rappis2,unpdel,y,yint,yj,yo,yo1,yppp,ypr |
|---|
| 44 | REAL yprimin,pisnmax,ydeg |
|---|
| 45 | INTEGER i,j,it,iter |
|---|
| 46 | c |
|---|
| 47 | c |
|---|
| 48 | c ------------------------------------------------------------------------- |
|---|
| 49 | |
|---|
| 50 | c On veut calculer y(Y) et y'(Y) , avec Y ,latitude de travail =f(j), |
|---|
| 51 | c et y , latitude vraie . |
|---|
| 52 | c |
|---|
| 53 | c Pour cela ,on se donne pour Y'(y) une fonction hyperbolique calculee |
|---|
| 54 | c sur ( 1+ nmax ) points ( plus nmax est grand, plus la precision pour |
|---|
| 55 | c y(Y) sera grande ) , on integre Y'(y) sur ces ( 1+nmax) pts pour |
|---|
| 56 | c obtenir Y(y) apres normalisation , puis on aura y(Y) en resolvant |
|---|
| 57 | c l'equation Y(yj) = Yj par Newton_Raphson , avec Yj = f(j) = |
|---|
| 58 | c pi/2 + pi/jjm ( 1.- j) pour les scalaires ( j =1, jjm+1 ) |
|---|
| 59 | c On en tire donc yj qui correspond a la ligne j . |
|---|
| 60 | c On a ensuite pour les derivees : y'(Y) = 1./ Y'(y) * ( pi/ jjm ) |
|---|
| 61 | c |
|---|
| 62 | c On calculera y(Y) et y'(Y) successivement aux latitudes de U |
|---|
| 63 | c , V , U + 0.25 , U + 0.75 . |
|---|
| 64 | c |
|---|
| 65 | c Notations ; |
|---|
| 66 | c ------------- |
|---|
| 67 | c yi represente les nmax+1 latitudes de depart ( -pi/2 a pi/2 ) |
|---|
| 68 | c yyy(nmax) repres. l'integrale de Y'(y) de - pi/2 a pi/2 |
|---|
| 69 | c yprim represente Y'(y) |
|---|
| 70 | c yf represente Y (y) |
|---|
| 71 | c |
|---|
| 72 | c yo represente la latitude du centre du zoom |
|---|
| 73 | c delta represente la demi largeur du zoom |
|---|
| 74 | c |
|---|
| 75 | c |
|---|
| 76 | c N.B : La fonction hyperbolique Y'(y) est definie comme suit : |
|---|
| 77 | c ------ |
|---|
| 78 | c |
|---|
| 79 | c Y' = TANH ( ( psio - psi )/psi*(1- psi) ) avec : |
|---|
| 80 | c |
|---|
| 81 | c psi = ( y - yo ) / ( pi/2 - yo ) si yo < y < pi/2 |
|---|
| 82 | c psi = (yo - y )/ ( pi/2 + yo ) si - pi/2 < y < yo |
|---|
| 83 | c |
|---|
| 84 | c |
|---|
| 85 | c .... yo represente la latitude du centre du zoom ..... |
|---|
| 86 | c |
|---|
| 87 | c .... psi0 et delta sont des parametres variant entre 0 et 1 .... |
|---|
| 88 | c---------------------------------------------------------------------- |
|---|
| 89 | c |
|---|
| 90 | cc psi0 = 0.3 et delta = 0.5 sont les valeurs qui conviennent |
|---|
| 91 | c pour *** iim = 96 et jjm = 72 *** , d'apres des tests sur |
|---|
| 92 | c les valeurs du Laplacien itere ( 2 fois) , analytiques et cal- |
|---|
| 93 | c culees. Pour d'autres valeurs de iim et jjm , on peut tester |
|---|
| 94 | c les nouvelles valeurs de psi0 et delta en lancant testfxyhyp . |
|---|
| 95 | c |
|---|
| 96 | c ---------------------------------------------------------------------- |
|---|
| 97 | c Fxyhp est appele par inigeom a la place de fxynew ( fonction |
|---|
| 98 | c a derivee sinusoidale ) si la variable logique fxyhypb , lue par |
|---|
| 99 | c defrun_new est . TRUE. . |
|---|
| 100 | c ---------------------------------------------------------------------- |
|---|
| 101 | c |
|---|
| 102 | |
|---|
| 103 | pi = 2.*ASIN( 1. ) |
|---|
| 104 | pis2 = pi/2. |
|---|
| 105 | yo = clat * pi/180. |
|---|
| 106 | psi0 = 0.3 |
|---|
| 107 | delta = 0.5 |
|---|
| 108 | unpdel = 1. + delta |
|---|
| 109 | eps = .1e-7 |
|---|
| 110 | c |
|---|
| 111 | PRINT *,' ----------------------------------------------------' |
|---|
| 112 | PRINT 1, psi0,delta,clat |
|---|
| 113 | PRINT *,' ----------------------------------------------------' |
|---|
| 114 | |
|---|
| 115 | pisnmax = pi/FLOAT(nmax) |
|---|
| 116 | yi(0) = - pis2 |
|---|
| 117 | |
|---|
| 118 | DO i = 1, nmax |
|---|
| 119 | yi(i)= - pis2 + FLOAT(i)*pisnmax |
|---|
| 120 | ENDDO |
|---|
| 121 | |
|---|
| 122 | yyy(0)= 0. |
|---|
| 123 | |
|---|
| 124 | c |
|---|
| 125 | c *************************************************************** |
|---|
| 126 | c ************ Cas ou yo < 0. ******************** |
|---|
| 127 | c *************************************************************** |
|---|
| 128 | c |
|---|
| 129 | IF( yo.LT. 0. ) THEN |
|---|
| 130 | c |
|---|
| 131 | pis2pyo = - pis2 + yo |
|---|
| 132 | rappis2 = ( -pis2-yo)/pis2pyo |
|---|
| 133 | rappis2 = rappis2 * rappis2 |
|---|
| 134 | c |
|---|
| 135 | DO 5 i = 1, nmax |
|---|
| 136 | |
|---|
| 137 | y = 0.5 * ( yi(i)+ yi(i-1) ) |
|---|
| 138 | IF( y.LT.- pis2.OR.y.GT.pis2 ) STOP 0 |
|---|
| 139 | |
|---|
| 140 | IF(y.GT.yo) THEN |
|---|
| 141 | psi = (y-yo)/(pis2-yo) |
|---|
| 142 | ypr = TANH( (psi0 -psi)/(psi*(1.-psi)) ) |
|---|
| 143 | |
|---|
| 144 | ELSE |
|---|
| 145 | psi = (yo-y)/(pis2+yo) |
|---|
| 146 | ypr = 1.-rappis2 + rappis2* TANH ( (psi0 -psi)/(psi*(1.-psi)) ) |
|---|
| 147 | ENDIF |
|---|
| 148 | |
|---|
| 149 | yyy(i)= yyy(i-1)+ ypr*(yi(i)-yi(i-1)) |
|---|
| 150 | |
|---|
| 151 | 5 CONTINUE |
|---|
| 152 | |
|---|
| 153 | den = unpdel + yyy(nmax)/pi |
|---|
| 154 | |
|---|
| 155 | yyy(0)= 0. |
|---|
| 156 | yf(0) = - pis2 |
|---|
| 157 | |
|---|
| 158 | yprim(nmax) = -1. |
|---|
| 159 | yprim(0) = 1.- rappis2 - rappis2 |
|---|
| 160 | |
|---|
| 161 | DO 12 i = 0 ,nmax |
|---|
| 162 | y = yi(i) |
|---|
| 163 | IF(y.LT.-pis2.OR.y.GT.pis2) STOP 1 |
|---|
| 164 | |
|---|
| 165 | IF(y.LT.yo) THEN |
|---|
| 166 | psi = (yo -y)/(pis2+yo) |
|---|
| 167 | IF(i.NE.0.AND.i.NE.nmax) THEN |
|---|
| 168 | yprim(i) = 1.- rappis2 + rappis2 |
|---|
| 169 | * * TANH( (psi0 -psi)/(psi*(1.-psi)) ) |
|---|
| 170 | ENDIF |
|---|
| 171 | ELSE |
|---|
| 172 | psi = (y-yo)/(pis2-yo) |
|---|
| 173 | IF(i.NE.0.AND.i.NE.nmax) THEN |
|---|
| 174 | yprim(i) = TANH( (psi0 -psi)/(psi*(1.-psi)) ) |
|---|
| 175 | ENDIF |
|---|
| 176 | ENDIF |
|---|
| 177 | |
|---|
| 178 | c |
|---|
| 179 | c ..................................................... |
|---|
| 180 | c ....... Normalisation de Y'(y) et de Y(y) ..... |
|---|
| 181 | c ..................................................... |
|---|
| 182 | c |
|---|
| 183 | yprim(i)= ( yprim(i)+ unpdel )/ den |
|---|
| 184 | yf(i)= (yyy(i)+ unpdel*( yi(i)+ pis2))/ den - pis2 |
|---|
| 185 | c |
|---|
| 186 | 12 CONTINUE |
|---|
| 187 | c |
|---|
| 188 | ELSE |
|---|
| 189 | c |
|---|
| 190 | c ************************************************************** |
|---|
| 191 | c ********** Cas ou yo >= 0 . ************** |
|---|
| 192 | c ************************************************************** |
|---|
| 193 | c |
|---|
| 194 | pis2pyo = pis2 + yo |
|---|
| 195 | rappis2 = ( pis2 - yo )/pis2pyo |
|---|
| 196 | rappis2 = rappis2 * rappis2 |
|---|
| 197 | c |
|---|
| 198 | |
|---|
| 199 | DO 13 i = 1, nmax |
|---|
| 200 | y= 0.5 * ( yi(i)+ yi(i-1) ) |
|---|
| 201 | IF(y.LT.- pis2.OR.y.GT.pis2) STOP 2 |
|---|
| 202 | |
|---|
| 203 | IF(y.LT.yo) THEN |
|---|
| 204 | psi = (yo -y)/(pis2+yo) |
|---|
| 205 | ypr = TANH( (psi0 - psi)/(psi*(1.-psi)) ) |
|---|
| 206 | |
|---|
| 207 | ELSE |
|---|
| 208 | psi = (y- yo)/(pis2-yo) |
|---|
| 209 | ypr = 1.-rappis2 + rappis2* TANH ( (psi0 -psi)/(psi*(1.-psi)) ) |
|---|
| 210 | |
|---|
| 211 | ENDIF |
|---|
| 212 | |
|---|
| 213 | yyy(i)= yyy(i-1)+ ypr*(yi(i)-yi(i-1)) |
|---|
| 214 | |
|---|
| 215 | 13 CONTINUE |
|---|
| 216 | |
|---|
| 217 | den = unpdel + yyy(nmax)/pi |
|---|
| 218 | |
|---|
| 219 | yyy(0)= 0. |
|---|
| 220 | yf(0) = -pis2 |
|---|
| 221 | |
|---|
| 222 | yprim(0) = -1. |
|---|
| 223 | yprim(nmax)= 1.- rappis2 - rappis2 |
|---|
| 224 | |
|---|
| 225 | DO 15 i = 0,nmax |
|---|
| 226 | y = yi(i) |
|---|
| 227 | IF(y.LT.-pis2.OR.y.GT.pis2) STOP 3 |
|---|
| 228 | |
|---|
| 229 | IF(y.LT.yo) THEN |
|---|
| 230 | psi = (yo -y)/(pis2+yo) |
|---|
| 231 | IF(i.NE.0.AND.i.NE.nmax) |
|---|
| 232 | * yprim(i) = TANH( (psi0 -psi)/(psi*(1.-psi)) ) |
|---|
| 233 | ELSE |
|---|
| 234 | psi = (y- yo)/(pis2-yo) |
|---|
| 235 | IF(i.NE.0.AND.i.NE.nmax) |
|---|
| 236 | * yprim(i) = 1.- rappis2 + rappis2 |
|---|
| 237 | * * TANH( (psi0 -psi)/(psi*(1.-psi)) ) |
|---|
| 238 | ENDIF |
|---|
| 239 | c ..................................................... |
|---|
| 240 | c |
|---|
| 241 | c ....... Normalisation de Y'(y) et de Y(y) ..... |
|---|
| 242 | c ..................................................... |
|---|
| 243 | c |
|---|
| 244 | yprim(i)= ( yprim(i)+ unpdel )/ den |
|---|
| 245 | yf(i)= (yyy(i)+ unpdel*( yi(i)+ pis2))/ den - pis2 |
|---|
| 246 | |
|---|
| 247 | 15 CONTINUE |
|---|
| 248 | |
|---|
| 249 | ENDIF |
|---|
| 250 | c |
|---|
| 251 | c *************************************************************** |
|---|
| 252 | |
|---|
| 253 | yf(0) = - pis2 |
|---|
| 254 | yf(nmax)= pis2 |
|---|
| 255 | |
|---|
| 256 | pijjm = pi/FLOAT(jjm) |
|---|
| 257 | |
|---|
| 258 | c ............................................................... |
|---|
| 259 | c |
|---|
| 260 | c Calculs de y(Y) et y'(Y) pour les latitudes de U ..... |
|---|
| 261 | c |
|---|
| 262 | PRINT *,' ***************************** ' |
|---|
| 263 | PRINT *,' *** Calcul de rlatu ***** ' |
|---|
| 264 | PRINT *,' ***************************** ' |
|---|
| 265 | c |
|---|
| 266 | c ............................................................... |
|---|
| 267 | c |
|---|
| 268 | PRINT *,' *** j rlatu yprimu **** ' |
|---|
| 269 | c |
|---|
| 270 | DO 1500 j=1, jjm +1 |
|---|
| 271 | |
|---|
| 272 | phij = pis2 + pijjm* ( 1.-Float(j) ) |
|---|
| 273 | yo1 = 0. |
|---|
| 274 | yj = yo1 |
|---|
| 275 | c |
|---|
| 276 | DO 500 iter = 1,300 |
|---|
| 277 | |
|---|
| 278 | DO 250 it = nmax,0,-1 |
|---|
| 279 | IF( yj.GE.yi(it)) GO TO 350 |
|---|
| 280 | 250 CONTINUE |
|---|
| 281 | |
|---|
| 282 | it= 0 |
|---|
| 283 | yj = yi(it) |
|---|
| 284 | |
|---|
| 285 | 350 CONTINUE |
|---|
| 286 | IF(it.EQ.nmax) THEN |
|---|
| 287 | yint = yf(nmax) |
|---|
| 288 | yprimin = yprim(nmax) |
|---|
| 289 | ELSE |
|---|
| 290 | c ................................................................. |
|---|
| 291 | c .... Interpolation entre yi(it) et yi(it+1) pour avoir Y(yj) |
|---|
| 292 | c ..... et Y'(yj) ..... |
|---|
| 293 | c ................................................................. |
|---|
| 294 | yint = (yf(it+1)-yf(it))/(yi(it+1)-yi(it))* ( yj-yi(it) ) + |
|---|
| 295 | + yf(it) |
|---|
| 296 | yprimin=(yprim(it+1)-yprim(it))/(yi(it+1)-yi(it))* ( yj-yi(it) ) + |
|---|
| 297 | + yprim(it) |
|---|
| 298 | ENDIF |
|---|
| 299 | yj = yj - (yint-phij)/yprimin |
|---|
| 300 | |
|---|
| 301 | IF( ABS(yj-yo1).LE.eps) GO TO 550 |
|---|
| 302 | yo1 = yj |
|---|
| 303 | c |
|---|
| 304 | 500 CONTINUE |
|---|
| 305 | PRINT *,' *** PAS DE SOLUTION **** ',j,phij |
|---|
| 306 | STOP 4 |
|---|
| 307 | 550 CONTINUE |
|---|
| 308 | ydeg = yj * 180./pi |
|---|
| 309 | rlatu(j) = yj |
|---|
| 310 | IF(j.EQ.1) rlatu(j) = pis2 |
|---|
| 311 | c |
|---|
| 312 | c ..... calcul de yprimu ..... |
|---|
| 313 | c |
|---|
| 314 | |
|---|
| 315 | IF(j.EQ.1.OR.j.EQ.jjm +1) GO TO 1500 |
|---|
| 316 | IF(yj.LT.yo) THEN |
|---|
| 317 | psi = (yo -yj)/(pis2+yo) |
|---|
| 318 | ELSE |
|---|
| 319 | psi = (yj-yo)/(pis2-yo) |
|---|
| 320 | ENDIF |
|---|
| 321 | |
|---|
| 322 | DO it = nmax,0,-1 |
|---|
| 323 | IF( yj.GE.yi(it)) GO TO 555 |
|---|
| 324 | ENDDO |
|---|
| 325 | 555 CONTINUE |
|---|
| 326 | IF(it.EQ.nmax) THEN |
|---|
| 327 | PRINT *,' IT = nmax Pbs ..',it,j |
|---|
| 328 | STOP 5 |
|---|
| 329 | ENDIF |
|---|
| 330 | yppp=(yprim(it+1)-yprim(it))/(yi(it+1)-yi(it))*(yj-yi(it)) + |
|---|
| 331 | * yprim(it) |
|---|
| 332 | yprimu(j) = pijjm/yppp |
|---|
| 333 | c |
|---|
| 334 | PRINT 3000,j,ydeg,yprimu(j) |
|---|
| 335 | 1500 CONTINUE |
|---|
| 336 | rlatu(1) = pis2 |
|---|
| 337 | rlatu(jjm+1) = - pis2 |
|---|
| 338 | c |
|---|
| 339 | c ............................................. |
|---|
| 340 | c |
|---|
| 341 | c ..... Calculs pour rlatv et yprimv ..... |
|---|
| 342 | c ............................................. |
|---|
| 343 | c |
|---|
| 344 | |
|---|
| 345 | PRINT *,' ***************************** ' |
|---|
| 346 | PRINT *,' **** Calcul de rlatv ***** ' |
|---|
| 347 | PRINT *,' ***************************** ' |
|---|
| 348 | PRINT *,' *** j rlatv yprimv **** ' |
|---|
| 349 | |
|---|
| 350 | DO 1600 j=1, jjm |
|---|
| 351 | |
|---|
| 352 | phij = pis2+ pijjm*(1.-(FLOAT(j)+0.5)) |
|---|
| 353 | yo1 = 0. |
|---|
| 354 | yj = yo1 |
|---|
| 355 | |
|---|
| 356 | DO 1550 iter =1,300 |
|---|
| 357 | |
|---|
| 358 | DO 1520 it = nmax,0,-1 |
|---|
| 359 | IF( yj.GE.yi(it)) GO TO 1530 |
|---|
| 360 | 1520 CONTINUE |
|---|
| 361 | c |
|---|
| 362 | it = 0 |
|---|
| 363 | yj = yi(it) |
|---|
| 364 | 1530 CONTINUE |
|---|
| 365 | IF(it.EQ.nmax) THEN |
|---|
| 366 | yint = yf(nmax) |
|---|
| 367 | yprimin = yprim(nmax) |
|---|
| 368 | ELSE |
|---|
| 369 | yint= (yf(it+1)-yf(it))/(yi(it+1)-yi(it)) * ( yj-yi(it) ) + |
|---|
| 370 | + yf(it) |
|---|
| 371 | yprimin=(yprim(it+1)-yprim(it))/(yi(it+1)-yi(it))* ( yj-yi(it) ) + |
|---|
| 372 | + yprim(it) |
|---|
| 373 | ENDIF |
|---|
| 374 | yj = yj - (yint-phij)/yprimin |
|---|
| 375 | |
|---|
| 376 | IF( ABS(yj-yo1).LE.eps ) GO TO 1560 |
|---|
| 377 | yo1 = yj |
|---|
| 378 | c |
|---|
| 379 | 1550 CONTINUE |
|---|
| 380 | PRINT *,' *** PAS DE SOLUTION *** ',j,phij |
|---|
| 381 | 1560 CONTINUE |
|---|
| 382 | ydeg = yj * 180./pi |
|---|
| 383 | rlatv(j)= yj |
|---|
| 384 | c |
|---|
| 385 | c ..... calcul de yprimv ..... |
|---|
| 386 | c |
|---|
| 387 | IF(yj.LT.-pis2.OR.yj.GT.pis2) STOP 6 |
|---|
| 388 | |
|---|
| 389 | IF(yj.LT.yo) THEN |
|---|
| 390 | psi = (yo -yj)/(pis2+yo) |
|---|
| 391 | ELSE |
|---|
| 392 | psi = (yj-yo)/(pis2-yo) |
|---|
| 393 | ENDIF |
|---|
| 394 | DO it = nmax,0,-1 |
|---|
| 395 | IF( yj.GE.yi(it)) GO TO 1570 |
|---|
| 396 | ENDDO |
|---|
| 397 | 1570 CONTINUE |
|---|
| 398 | IF(it.EQ.nmax) THEN |
|---|
| 399 | PRINT *,' IT = nmax Pbs ..',it,j |
|---|
| 400 | STOP 7 |
|---|
| 401 | ENDIF |
|---|
| 402 | yppp = (yprim(it+1)-yprim(it))/(yi(it+1)-yi(it))*(yj-yi(it)) + |
|---|
| 403 | * yprim(it) |
|---|
| 404 | yprimv(j)= pijjm/ yppp |
|---|
| 405 | c |
|---|
| 406 | PRINT 3000,j,ydeg,yprimv(j) |
|---|
| 407 | |
|---|
| 408 | 1600 CONTINUE |
|---|
| 409 | |
|---|
| 410 | c |
|---|
| 411 | c .................................................. |
|---|
| 412 | c |
|---|
| 413 | PRINT *,' ************************************ ' |
|---|
| 414 | PRINT *,' **** Calcul pour f( j + 0.25 ) **** ' |
|---|
| 415 | PRINT *,' ************************************ ' |
|---|
| 416 | PRINT *,' *** j rlatu1 yprimu1 **** ' |
|---|
| 417 | c |
|---|
| 418 | c .................................................. |
|---|
| 419 | c |
|---|
| 420 | DO 1800 j=1, jjm |
|---|
| 421 | |
|---|
| 422 | phij = pis2+ pijjm*( 1.- (FLOAT(j) + 0.25) ) |
|---|
| 423 | yo1 = 0. |
|---|
| 424 | yj = yo1 |
|---|
| 425 | C |
|---|
| 426 | |
|---|
| 427 | DO 1620 iter = 1, 500 |
|---|
| 428 | |
|---|
| 429 | DO 1610 it = nmax,0,-1 |
|---|
| 430 | IF( yj.GE.yi(it)) GO TO 1615 |
|---|
| 431 | 1610 CONTINUE |
|---|
| 432 | c |
|---|
| 433 | it = 0 |
|---|
| 434 | yj = yi(it) |
|---|
| 435 | 1615 CONTINUE |
|---|
| 436 | IF(it.EQ.nmax) THEN |
|---|
| 437 | yint = yf(nmax) |
|---|
| 438 | yprimin = yprim(nmax) |
|---|
| 439 | ELSE |
|---|
| 440 | yint= (yf(it+1)-yf(it))/(yi(it+1)-yi(it))*(yj-yi(it))+ |
|---|
| 441 | * yf(it) |
|---|
| 442 | yprimin=(yprim(it+1)-yprim(it))/(yi(it+1)-yi(it))*(yj-yi(it))+ |
|---|
| 443 | * yprim(it) |
|---|
| 444 | ENDIF |
|---|
| 445 | yj = yj - (yint-phij)/yprimin |
|---|
| 446 | c IF( (- pis2-yj)*(yj- pis2 ).LT.0.)PRINT *,' Newton hors limites' |
|---|
| 447 | |
|---|
| 448 | IF( ABS(yj-yo1).LE.eps ) GO TO 1625 |
|---|
| 449 | yo1 = yj |
|---|
| 450 | c |
|---|
| 451 | 1620 CONTINUE |
|---|
| 452 | PRINT *,' *** PAS DE SOLUTION **** ',j,phij |
|---|
| 453 | 1625 CONTINUE |
|---|
| 454 | ydeg = yj * 180./pi |
|---|
| 455 | rlatu1(j)= yj |
|---|
| 456 | c |
|---|
| 457 | c ..... calcul de yprimu1 ..... |
|---|
| 458 | c |
|---|
| 459 | IF( yj.LT.-pis2.OR.yj.GT.pis2 ) STOP 8 |
|---|
| 460 | c |
|---|
| 461 | IF(yj.LT.yo) THEN |
|---|
| 462 | psi = (yo -yj)/(pis2+yo) |
|---|
| 463 | ELSE |
|---|
| 464 | psi = (yj-yo)/(pis2-yo) |
|---|
| 465 | ENDIF |
|---|
| 466 | DO it = nmax,0,-1 |
|---|
| 467 | IF( yj.GE.yi(it)) GO TO 1630 |
|---|
| 468 | ENDDO |
|---|
| 469 | 1630 CONTINUE |
|---|
| 470 | IF(it.EQ.nmax) THEN |
|---|
| 471 | PRINT *,' IT = nmax Pbs ..',it,j |
|---|
| 472 | STOP 9 |
|---|
| 473 | ENDIF |
|---|
| 474 | yppp = (yprim(it+1)-yprim(it))/(yi(it+1)-yi(it))*(yj-yi(it)) + |
|---|
| 475 | * yprim(it) |
|---|
| 476 | yprimu1(j) = pijjm/yppp |
|---|
| 477 | c |
|---|
| 478 | PRINT 3000,j,ydeg,yprimu1(j) |
|---|
| 479 | |
|---|
| 480 | 1800 CONTINUE |
|---|
| 481 | |
|---|
| 482 | c .............................................. |
|---|
| 483 | c |
|---|
| 484 | PRINT *,' ********************************** ' |
|---|
| 485 | PRINT *,' **** Calcul pour f( j + 3/4 ) *** ' |
|---|
| 486 | PRINT *,' ********************************** ' |
|---|
| 487 | PRINT *,' *** j rlatu2 yprimu2 *** ' |
|---|
| 488 | c .............................................. |
|---|
| 489 | c |
|---|
| 490 | DO 1900 j=1, jjm |
|---|
| 491 | |
|---|
| 492 | phij = pis2+ pijjm*(1.-(FLOAT(j)+0.75)) |
|---|
| 493 | yo1 = 0. |
|---|
| 494 | yj = yo1 |
|---|
| 495 | |
|---|
| 496 | DO 1700 iter =1,500 |
|---|
| 497 | |
|---|
| 498 | DO 1680 it = nmax,0,-1 |
|---|
| 499 | IF( yj.GE.yi(it)) GO TO 1690 |
|---|
| 500 | 1680 CONTINUE |
|---|
| 501 | it = 0 |
|---|
| 502 | yj = yi(it) |
|---|
| 503 | 1690 CONTINUE |
|---|
| 504 | IF(it.EQ.nmax) THEN |
|---|
| 505 | yint = yf(nmax) |
|---|
| 506 | yprimin = yprim(nmax) |
|---|
| 507 | ELSE |
|---|
| 508 | yint= (yf(it+1)-yf(it))/(yi(it+1)-yi(it)) * ( yj -yi(it) ) + |
|---|
| 509 | + yf(it) |
|---|
| 510 | yprimin=(yprim(it+1)-yprim(it))/(yi(it+1)-yi(it))* ( yj-yi(it) ) + |
|---|
| 511 | + yprim(it) |
|---|
| 512 | ENDIF |
|---|
| 513 | yj = yj - (yint-phij)/yprimin |
|---|
| 514 | |
|---|
| 515 | IF( ABS(yj-yo1).LE.eps ) GO TO 1705 |
|---|
| 516 | yo1 = yj |
|---|
| 517 | c |
|---|
| 518 | 1700 CONTINUE |
|---|
| 519 | PRINT *,' PAS DE SOLUTION ',j,phij |
|---|
| 520 | STOP 10 |
|---|
| 521 | 1705 CONTINUE |
|---|
| 522 | ydeg = yj * 180./pi |
|---|
| 523 | rlatu2(j) = yj |
|---|
| 524 | c |
|---|
| 525 | c ..... calcul de yprimu2 ..... |
|---|
| 526 | c |
|---|
| 527 | IF(yj.LT.- pis2.OR.yj.GT.pis2) STOP 11 |
|---|
| 528 | |
|---|
| 529 | IF(yj.LT.yo) THEN |
|---|
| 530 | psi = (yo -yj)/(pis2+yo) |
|---|
| 531 | ELSE |
|---|
| 532 | psi = (yj-yo)/(pis2-yo) |
|---|
| 533 | ENDIF |
|---|
| 534 | c |
|---|
| 535 | DO it = nmax,0,-1 |
|---|
| 536 | IF( yj.GE.yi(it)) GO TO 1710 |
|---|
| 537 | ENDDO |
|---|
| 538 | 1710 CONTINUE |
|---|
| 539 | IF(it.EQ.nmax) THEN |
|---|
| 540 | PRINT *,' IT = nmax Pbs ..',it,j |
|---|
| 541 | STOP 12 |
|---|
| 542 | ENDIF |
|---|
| 543 | yppp =(yprim(it+1)-yprim(it))/(yi(it+1)-yi(it))* ( yj - yi(it) ) + |
|---|
| 544 | + yprim(it) |
|---|
| 545 | yprimu2(j) = pijjm/yppp |
|---|
| 546 | c |
|---|
| 547 | PRINT 3000,j,ydeg,yprimu2(j) |
|---|
| 548 | |
|---|
| 549 | 1900 CONTINUE |
|---|
| 550 | c |
|---|
| 551 | PRINT *,' *** TEST DE COHERENCE DANS FXYP **** ' |
|---|
| 552 | c |
|---|
| 553 | DO j = 1, jjm |
|---|
| 554 | c |
|---|
| 555 | IF(rlatu1(j).LE.rlatu2(j)) THEN |
|---|
| 556 | PRINT *,' Attention ! rlatu1 < rlatu2 ',rlatu1(j), rlatu2(j),j |
|---|
| 557 | STOP 13 |
|---|
| 558 | ENDIF |
|---|
| 559 | c |
|---|
| 560 | IF(rlatu2(j).LE.rlatu(j+1)) THEN |
|---|
| 561 | PRINT *,' Attention ! rlatu2 < rlatup1 ',rlatu2(j),rlatu(j+1),j |
|---|
| 562 | STOP 14 |
|---|
| 563 | ENDIF |
|---|
| 564 | c |
|---|
| 565 | IF(rlatu(j).LE.rlatu1(j)) THEN |
|---|
| 566 | PRINT *,' rlatu < rlatu1 ',rlatu(j),rlatu1(j),j |
|---|
| 567 | STOP 15 |
|---|
| 568 | ENDIF |
|---|
| 569 | c |
|---|
| 570 | IF(rlatv(j).LE.rlatu2(j)) THEN |
|---|
| 571 | PRINT *,' rlatv > rlatu2 ',rlatv(j),rlatu2(j),j |
|---|
| 572 | STOP 16 |
|---|
| 573 | ENDIF |
|---|
| 574 | c |
|---|
| 575 | IF(rlatv(j).ge.rlatu1(j)) THEN |
|---|
| 576 | PRINT *,' rlatv < rlatu1 ',rlatv(j),rlatu1(j),j |
|---|
| 577 | STOP 17 |
|---|
| 578 | ENDIF |
|---|
| 579 | c |
|---|
| 580 | IF(rlatv(j).ge.rlatu(j)) THEN |
|---|
| 581 | PRINT *,' rlatv > rlatu ',rlatv(j),rlatu(j),j |
|---|
| 582 | STOP 18 |
|---|
| 583 | ENDIF |
|---|
| 584 | c |
|---|
| 585 | ENDDO |
|---|
| 586 | c |
|---|
| 587 | PRINT *,' *** Les latitudes yprimu,yprimv,yprimu1,yprimu2 ', |
|---|
| 588 | * ' sont coherentes entre elles ! *** ' |
|---|
| 589 | |
|---|
| 590 | c |
|---|
| 591 | 1 FORMAT( 3x,'Fxyhyp. psi0 delta clat ',3f7.2) |
|---|
| 592 | 3000 FORMAT(2x,i4,2x,f8.2,f8.3) |
|---|
| 593 | RETURN |
|---|
| 594 | END |
|---|
