Changeset 203 for LMDZ.3.3/trunk/libf
- Timestamp:
- Apr 12, 2001, 11:30:31 AM (24 years ago)
- Location:
- LMDZ.3.3/trunk/libf/dyn3d
- Files:
-
- 1 added
- 8 edited
-
defrun_new.F (modified) (18 diffs)
-
dynetat0.F (modified) (2 diffs)
-
dynredem.F (modified) (3 diffs)
-
fxhyp.F (modified) (10 diffs)
-
fxyhyper.F (modified) (7 diffs)
-
fyhyp.F (modified) (11 diffs)
-
heavyside.F (added)
-
inigeom.F (modified) (13 diffs)
-
serre.h (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ.3.3/trunk/libf/dyn3d/defrun_new.F
r2 r203 1 c 2 c $Header 3 c 1 4 SUBROUTINE defrun_new( tapedef, etatinit, clesphy0 ) 2 5 c … … 36 39 INTEGER tapeout 37 40 REAL clonn,clatt,grossismxx,grossismyy 38 REAL dzoomxx,dzoomyy 41 REAL dzoomxx,dzoomyy,tauxx,tauyy 39 42 LOGICAL fxyhypbb, ysinuss 40 43 INTEGER i … … 258 261 WRITE(tapeout,*) clonn 259 262 IF( ABS(clon - clonn).GE. 0.001 ) THEN 260 PRINT *,' La valeur de clon passee par run.def est differente de261 * celle lue sur le fichier start '263 WRITE(tapeout,*) ' La valeur de clon passee par run.def est diffe 264 *rente de celle lue sur le fichier start ' 262 265 STOP 263 266 ENDIF 264 c265 267 c 266 268 READ (tapedef,9001) ch1,ch4 … … 270 272 271 273 IF( ABS(clat - clatt).GE. 0.001 ) THEN 272 PRINT *,' La valeur de clat passee par run.def est differente de273 * celle lue sur le fichier start '274 WRITE(tapeout,*) ' La valeur de clat passee par run.def est diffe 275 *rente de celle lue sur le fichier start ' 274 276 STOP 275 277 ENDIF … … 281 283 282 284 IF( ABS(grossismx - grossismxx).GE. 0.001 ) THEN 283 PRINT *,' La valeur de grossismx passee par run.def est differente284 *de celle lue sur le fichier start '285 WRITE(tapeout,*) ' La valeur de grossismx passee par run.def est 286 , differente de celle lue sur le fichier start ' 285 287 STOP 286 288 ENDIF … … 292 294 293 295 IF( ABS(grossismy - grossismyy).GE. 0.001 ) THEN 294 PRINT *,' La valeur de grossismy passee par run.def est differen295 *te de celle lue sur le fichier start '296 WRITE(tapeout,*) ' La valeur de grossismy passee par run.def est 297 , differente de celle lue sur le fichier start ' 296 298 STOP 297 299 ENDIF 298 300 299 301 IF( grossismx.LT.1. ) THEN 300 PRINT *,' *** ATTENTION !! grossismx < 1 . *** '302 WRITE(tapeout,*) ' *** ATTENTION !! grossismx < 1 . *** ' 301 303 STOP 302 304 ELSE … … 306 308 307 309 IF( grossismy.LT.1. ) THEN 308 PRINT *,' *** ATTENTION !! grossismy < 1 . *** '310 WRITE(tapeout,*) ' *** ATTENTION !! grossismy < 1 . *** ' 309 311 STOP 310 312 ELSE … … 312 314 ENDIF 313 315 314 PRINT *,' alphax alphay defrun ',alphax,alphay315 316 c 316 317 c alphax et alphay sont les anciennes formulat. des grossissements … … 324 325 IF( .NOT.fxyhypb ) THEN 325 326 IF( fxyhypbb ) THEN 326 PRINT *,' ******** PBS DANS DEFRUN ******** '327 PRINT *,' *** fxyhypb lu sur le fichier start est F ',328 * 'alors qu il est T sur run.def ***'327 WRITE(tapeout,*) ' ******** PBS DANS DEFRUN ******** ' 328 WRITE(tapeout,*)' *** fxyhypb lu sur le fichier start est F' 329 *, ' alors qu il est T sur run.def ***' 329 330 STOP 330 331 ENDIF 331 332 ELSE 332 333 IF( .NOT.fxyhypbb ) THEN 333 PRINT *,' ******** PBS DANS DEFRUN ******** '334 PRINT *,' *** fxyhypb lu sur le fichier start est T ',335 * 'alors qu il est F sur run.def ****'334 WRITE(tapeout,*) ' ******** PBS DANS DEFRUN ******** ' 335 WRITE(tapeout,*)' *** fxyhypb lu sur le fichier start est t' 336 *, ' alors qu il est F sur run.def ***' 336 337 STOP 337 338 ENDIF … … 343 344 WRITE(tapeout,*) dzoomxx 344 345 345 IF( fxyhypb ) THEN346 IF( ABS(dzoomx - dzoomxx).GE. 0.001 ) THEN347 PRINT *,' La valeur de dzoomx passee par run.def est differente348 * de celle lue sur le fichier start '349 STOP350 ENDIF351 ENDIF352 353 346 READ (tapedef,9001) ch1,ch4 354 347 READ (tapedef,*) dzoomyy … … 356 349 WRITE(tapeout,*) dzoomyy 357 350 351 READ (tapedef,9001) ch1,ch4 352 READ (tapedef,*) tauxx 353 WRITE(tapeout,9001) ch1,'taux' 354 WRITE(tapeout,*) tauxx 355 356 READ (tapedef,9001) ch1,ch4 357 READ (tapedef,*) tauyy 358 WRITE(tapeout,9001) ch1,'tauy' 359 WRITE(tapeout,*) tauyy 360 358 361 IF( fxyhypb ) THEN 362 363 IF( ABS(dzoomx - dzoomxx).GE. 0.001 ) THEN 364 WRITE(tapeout,*)' La valeur de dzoomx passee par run.def est dif 365 *ferente de celle lue sur le fichier start ' 366 CALL ABORT 367 ENDIF 368 359 369 IF( ABS(dzoomy - dzoomyy).GE. 0.001 ) THEN 360 PRINT *,' La valeur de dzoomy passee par run.def est differente361 * de celle lue sur le fichier start '362 STOP370 WRITE(tapeout,*)' La valeur de dzoomy passee par run.def est dif 371 *ferente de celle lue sur le fichier start ' 372 CALL ABORT 363 373 ENDIF 374 375 IF( ABS(taux - tauxx).GE. 0.001 ) THEN 376 WRITE(6,*)' La valeur de taux passee par run.def est differente 377 * de celle lue sur le fichier start ' 378 CALL ABORT 379 ENDIF 380 381 IF( ABS(tauy - tauyy).GE. 0.001 ) THEN 382 WRITE(6,*)' La valeur de tauy passee par run.def est differente 383 * de celle lue sur le fichier start ' 384 CALL ABORT 385 ENDIF 386 364 387 ENDIF 365 388 … … 374 397 IF( .NOT.ysinus ) THEN 375 398 IF( ysinuss ) THEN 376 PRINT *,' ******** PBS DANS DEFRUN ******** '377 PRINT *,' *** ysinus lu sur le fichier start est F',378 * ' alors qu il est T sur run.def ***'399 WRITE(6,*) ' ******** PBS DANS DEFRUN ******** ' 400 WRITE(tapeout,*)'** ysinus lu sur le fichier start est F', 401 * ' alors qu il est T sur run.def ***' 379 402 STOP 380 403 ENDIF 381 404 ELSE 382 405 IF( .NOT.ysinuss ) THEN 383 PRINT *,' ******** PBS DANS DEFRUN ******** '384 PRINT *,' *** ysinus lu sur le fichier start est T',385 * 'alors qu il est F sur run.def ****'406 WRITE(6,*) ' ******** PBS DANS DEFRUN ******** ' 407 WRITE(tapeout,*)'** ysinus lu sur le fichier start est T', 408 * ' alors qu il est F sur run.def ***' 386 409 STOP 387 410 ENDIF … … 389 412 ENDIF 390 413 c 391 close(tapedef) 414 WRITE(6,*) ' alphax alphay defrun ',alphax,alphay 415 416 CLOSE(tapedef) 417 392 418 RETURN 393 419 c ............................................... … … 416 442 417 443 IF( grossismx.LT.1. ) THEN 418 PRINT *,'*** ATTENTION !! grossismx < 1 . *** '444 WRITE(tapeout,*) '*** ATTENTION !! grossismx < 1 . *** ' 419 445 STOP 420 446 ELSE … … 423 449 424 450 IF( grossismy.LT.1. ) THEN 425 PRINT *,' *** ATTENTION !! grossismy < 1 . *** '451 WRITE(tapeout,*) ' *** ATTENTION !! grossismy < 1 . *** ' 426 452 STOP 427 453 ELSE … … 429 455 ENDIF 430 456 431 PRINT *,' alphax alphay defrun ',alphax,alphay432 433 457 c 434 458 READ (tapedef,9001) ch1,ch4 … … 446 470 WRITE(tapeout,9001) ch1,'dzoomy' 447 471 WRITE(tapeout,*) dzoomy 448 c 472 473 READ (tapedef,9001) ch1,ch4 474 READ (tapedef,*) taux 475 WRITE(tapeout,9001) ch1,'taux' 476 WRITE(tapeout,*) taux 477 c 478 READ (tapedef,9001) ch1,ch4 479 READ (tapedef,*) tauy 480 WRITE(tapeout,9001) ch1,'tauy' 481 WRITE(tapeout,*) tauy 482 449 483 READ (tapedef,9001) ch1,ch4 450 484 READ (tapedef,*) ysinus … … 452 486 WRITE(tapeout,*) ysinus 453 487 488 WRITE(tapeout,*) ' alphax alphay defrun ',alphax,alphay 454 489 c 455 490 9000 FORMAT(3(/,a72)) 456 491 9001 FORMAT(/,a72,/,a12) 457 492 cc 458 close(tapedef) 493 CLOSE(tapedef) 494 459 495 RETURN 460 496 END -
LMDZ.3.3/trunk/libf/dyn3d/dynetat0.F
r2 r203 1 c 2 c $Header 3 c 1 4 SUBROUTINE dynetat0(fichnom,nq,vcov,ucov, 2 5 . teta,q,masse,ps,phis,time) … … 105 108 dzoomx = tab_cntrl(25) 106 109 dzoomy = tab_cntrl(26) 110 taux = tab_cntrl(28) 111 tauy = tab_cntrl(29) 107 112 ELSE 108 113 fxyhypb = . FALSE . -
LMDZ.3.3/trunk/libf/dyn3d/dynredem.F
r56 r203 1 c 2 c $Header 3 c 1 4 SUBROUTINE dynredem0(fichnom,idayref,anneeref,phis,nq) 2 5 USE IOIPSL … … 92 95 tab_cntrl(25) = dzoomx 93 96 tab_cntrl(26) = dzoomy 97 tab_cntrl(27) = 0. 98 tab_cntrl(28) = taux 99 tab_cntrl(29) = tauy 94 100 ELSE 95 101 tab_cntrl(24) = 0. … … 97 103 tab_cntrl(26) = dzoomy 98 104 tab_cntrl(27) = 0. 99 IF( ysinus ) tab_cntrl(27) = 1. 105 tab_cntrl(28) = 0. 106 tab_cntrl(29) = 0. 107 IF( ysinus ) tab_cntrl(27) = 1. 100 108 ENDIF 101 109 c -
LMDZ.3.3/trunk/libf/dyn3d/fxhyp.F
r2 r203 1 c 2 c $Header 3 c 1 4 SUBROUTINE fxhyp ( xzoomdeg,grossism,dzoom,tau , 2 5 , rlonm025,xprimm025,rlonv,xprimv,rlonu,xprimu,rlonp025,xprimp025) … … 11 14 c grossism etant le grossissement ( = 2 si 2 fois, = 3 si 3 fois,etc.) 12 15 c dzoom etant la distance totale de la zone du zoom 13 c tau la transition , normalement = 1 . 14 c 16 c tau la raideur de la transition de l'interieur a l'exterieur du zoom 17 c 18 c On doit avoir grossism x dzoom < pi ( radians ) , en longitude. 19 c ******************************************************************** 20 15 21 16 22 INTEGER nmax, nmax2 17 PARAMETER ( nmax = 50000, nmax2 = 2*nmax )23 PARAMETER ( nmax = 30000, nmax2 = 2*nmax ) 18 24 19 25 #include "dimensions.h" … … 23 29 c 24 30 REAL xzoomdeg,dzoom,tau,grossism 31 32 c ...... arguments de sortie ...... 33 25 34 REAL rlonm025(iip1),xprimm025(iip1),rlonv(iip1),xprimv(iip1), 26 35 , rlonu(iip1),xprimu(iip1),rlonp025(iip1),xprimp025(iip1) 27 36 28 c ...... arguments de sortie ...... 29 c 30 REAL xlon(iip1),xprimm(iip1),xuv 31 REAL xtild(0:nmax2) 32 REAL fhyp(0:nmax),ffdx(0:nmax),beta,Xprimt(0:nmax2) 33 REAL Xf(0:nmax2),xxpr(0:nmax2) 34 REAL xvrai(iip1),xxprim(iip1) 35 REAL pi,depi,epsilon,xzoom 37 c .... variables locales .... 38 c 39 REAL*8 xlon(iip1),xprimm(iip1),xuv 40 REAL*8 xtild(0:nmax2) 41 REAL*8 fhyp(0:nmax),ffdx,beta,Xprimt(0:nmax2) 42 REAL*8 Xf(0:nmax2),xxpr(0:nmax2) 43 REAL*8 xvrai(iip1),xxprim(iip1) 44 REAL*8 pi,depi,epsilon,xzoom,fa,fb 45 REAL*8 Xf1, Xfi , a0,a1,a2,a3,xi2 36 46 INTEGER i,it,ik,iter,ii,idif 37 REAL xi,xo1,xint,xmoy,xlon2,fxm,Xprimin38 REAL champmin,champmax47 REAL*8 xi,xo1,xmoy,xlon2,fxm,Xprimin 48 REAL*8 champmin,champmax 39 49 INTEGER is2 40 50 SAVE is2 41 REAL dlon1(iip1),dlon2(iip1),dlon3(iip1) 51 REAL*8 heavyside 52 EXTERNAL coefpoly,heavyside 42 53 43 54 pi = 2. * ASIN(1.) 44 55 depi = 2. * pi 45 epsilon = 1.e- 656 epsilon = 1.e-3 46 57 xzoom = xzoomdeg * pi/180. 47 58 48 59 WRITE(6,24) xzoomdeg,grossism,tau,dzoom 60 IF( dzoom.LT.1.) THEN 61 dzoom = dzoom * depi 62 ELSEIF( dzoom.LT. 25. ) THEN 63 WRITE(6,*) ' Le param. dzoomy pour fxhyp est trop petit ! L aug 64 ,menter et relancer ! ' 65 STOP 1 66 ELSE 67 dzoom = dzoom * pi/180. 68 ENDIF 49 69 50 70 DO i = 0, nmax2 51 xtild(i) = FLOAT(i) /nmax2 52 IF( xtild(i).EQ. 0.5 ) xtild(i) = xtild(i) + 1.e-6 53 ENDDO 54 55 DO i = 1, nmax 56 fhyp(i) = TANH ( ( xtild(i) - 0.5*(1.- dzoom) ) / 57 , ( tau * xtild(i) * ( 0.5 -xtild(i))) ) 58 ENDDO 59 60 fhyp( 0 ) = - 1. 61 fhyp( nmax ) = 1. 71 xtild(i) = - pi + FLOAT(i) * depi /nmax2 72 ENDDO 73 74 DO i = 0, nmax2 75 76 fa = tau* ( dzoom/2. - xtild(i) ) 77 fb = xtild(i) * ( pi - xtild(i) ) 78 79 80 IF( 200.* fb .LT. - fa ) THEN 81 fhyp ( i) = - 1. 82 ELSEIF( 200. * fb .LT. fa ) THEN 83 fhyp ( i) = 1. 84 ELSE 85 fhyp(i) = TANH ( fa/fb ) 86 ENDIF 87 88 IF ( xtild(i).EQ. 0. ) fhyp(i) = 1. 89 IF ( xtild(i).EQ. pi ) fhyp(i) = -1. 90 91 ENDDO 62 92 63 93 cc .... Calcul de beta .... 64 c 65 ffdx( 0 ) = 0. 66 67 DO i = 1, nmax 68 xmoy = 0.5 * ( xtild(i-1) + xtild( i ) ) 69 fxm = TANH ( ( xmoy - 0.5 * ( 1. - dzoom ) ) / 70 , ( tau * xmoy * ( 0.5 -xmoy)) ) 71 ffdx(i) = ffdx(i-1) + fxm * ( xtild(i) - xtild(i-1) ) 72 ENDDO 73 74 beta = ( grossism * ffdx(nmax) - 0.5 ) / ( ffdx(nmax) - 0.5 ) 94 c ............................ 95 96 ffdx = 0. 97 98 DO i = nmax +1,nmax2 99 100 xmoy = 0.5 * ( xtild(i-1) + xtild( i ) ) 101 fa = tau* ( dzoom/2. - xmoy ) 102 fb = xmoy * ( pi - xmoy ) 103 104 IF( 200.* fb .LT. - fa ) THEN 105 fxm = - 1. 106 ELSEIF( 200. * fb .LT. fa ) THEN 107 fxm = 1. 108 ELSE 109 fxm = TANH ( fa/fb ) 110 ENDIF 111 112 IF ( xmoy.EQ. 0. ) fhyp(i) = 1. 113 IF ( xmoy.EQ. pi ) fhyp(i) = -1. 114 ffdx = ffdx + fxm * ( xtild(i) - xtild(i-1) ) 115 116 ENDDO 117 118 beta = ( grossism * ffdx - pi ) / ( ffdx - pi ) 119 ccc WRITE(6,*) ' ** X beta **',beta,grossism 120 121 IF( 2.*beta - grossism.LE. 0.) THEN 122 WRITE(6,*) ' ** Attention ! La valeur beta calculee dans la rou 123 ,tine fxhyp est mauvaise ! ' 124 WRITE(6,*)'Modifier les valeurs de grossismx ,tau ou dzoomx ', 125 , ' et relancer ! *** ' 126 CALL ABORT 127 ENDIF 75 128 c 76 129 c ..... calcul de Xprimt ..... 77 130 c 78 131 79 DO i = 0, nmax132 DO i = nmax, nmax2 80 133 Xprimt(i) = beta + ( grossism - beta ) * fhyp(i) 81 134 ENDDO 82 135 c 83 DO i = 0, nmax136 DO i = nmax+1, nmax2 84 137 Xprimt( nmax2 - i ) = Xprimt( i ) 85 138 ENDDO … … 88 141 c ..... Calcul de Xf ........ 89 142 90 Xf(0) = 0. 91 DO i = 1, nmax 92 xmoy = 0.5 * ( xtild(i-1) + xtild( i ) ) 93 fxm = TANH ( ( xmoy - 0.5 * ( 1. - dzoom ) ) / 94 , ( tau * xmoy * ( 0.5 -xmoy)) ) 95 xxpr(i) = beta + ( grossism - beta ) * fxm 96 ENDDO 97 98 DO i = 1,nmax 99 xxpr(nmax2-i+1) = xxpr(i) 143 Xf(0) = - pi 144 145 DO i = nmax +1, nmax2 146 147 xmoy = 0.5 * ( xtild(i-1) + xtild( i ) ) 148 fa = tau* ( dzoom/2. - xmoy ) 149 fb = xmoy * ( pi - xmoy ) 150 151 IF( 200.* fb .LT. - fa ) THEN 152 fxm = - 1. 153 ELSEIF( 200. * fb .LT. fa ) THEN 154 fxm = 1. 155 ELSE 156 fxm = TANH ( fa/fb ) 157 ENDIF 158 159 IF ( xmoy.EQ. 0. ) fxm = 1. 160 IF ( xmoy.EQ. pi ) fxm = -1. 161 xxpr(i) = beta + ( grossism - beta ) * fxm 162 163 ENDDO 164 165 DO i = nmax+1, nmax2 166 xxpr(nmax2-i+1) = xxpr(i) 100 167 ENDDO 101 168 … … 103 170 Xf(i) = Xf(i-1) + xxpr(i) * ( xtild(i) - xtild(i-1) ) 104 171 ENDDO 105 do i=1,nmax2106 xf(i)=xf(i)/xf(nmax2)107 enddo108 109 110 PRINT *,' XF ',xf(0),xf(nmax),xf(nmax2)111 172 112 173 … … 130 191 ENDIF 131 192 193 xo1 = 0. 132 194 133 195 DO 1500 i = 1, iim 134 196 135 xlon2 = ( FLOAT(i) + xuv - 0.75) / FLOAT(iim) 136 137 xo1 = 0. 138 xi = xlon2 139 c 140 DO 500 iter = 1,300 141 197 xlon2 = - pi + (FLOAT(i) + xuv - 0.75) * depi / FLOAT(iim) 198 199 Xfi = xlon2 200 c 142 201 DO 250 it = nmax2,0,-1 143 IF( xi.GE.xtild(it)) GO TO 350202 IF( Xfi.GE.Xf(it)) GO TO 350 144 203 250 CONTINUE 145 204 146 205 it = 0 147 xi = xtild(it)148 206 149 207 350 CONTINUE 150 IF(it.EQ.nmax2) THEN 151 it = nmax2 -1 152 xf(it+1) = 1. 153 ENDIF 154 c ................................................................. 155 c .... Interpolation entre xi(it) et xi(it+1) pour avoir X(xi) 156 c ..... et X'(xi) ..... 157 c ................................................................. 158 159 xint = ( Xf(it+1)-Xf(it) ) / ( xtild(it+1)-xtild(it) ) * 160 + ( xi-xtild(it) ) + Xf(it) 161 Xprimin = ( Xprimt(it+1)-Xprimt(it) )/ ( xtild(it+1)-xtild(it) ) * 162 + ( xi-xtild(it) ) + Xprimt(it) 163 164 xi = xi - (xint-xlon2)/Xprimin 165 166 IF( ABS(xi-xo1).LE.epsilon) GO TO 550 167 xo1 = xi 168 c 208 209 c ...... Calcul de Xf(xi) ...... 210 c 211 xi = xtild(it) 212 213 IF(it.EQ.nmax2) THEN 214 it = nmax2 -1 215 Xf(it+1) = pi 216 ENDIF 217 c ..................................................................... 218 c 219 c Appel de la routine qui calcule les coefficients a0,a1,a2,a3 d'un 220 c polynome de degre 3 qui passe par les points (Xf(it),xtild(it) ) 221 c et (Xf(it+1),xtild(it+1) ) 222 223 CALL coefpoly ( Xf(it),Xf(it+1),Xprimt(it),Xprimt(it+1), 224 , xtild(it),xtild(it+1), a0, a1, a2, a3 ) 225 226 Xf1 = Xf(it) 227 Xprimin = a1 + 2.* a2 * xi + 3.*a3 * xi *xi 228 229 DO 500 iter = 1,300 230 xi = xi - ( Xf1 - Xfi )/ Xprimin 231 232 IF( ABS(xi-xo1).LE.epsilon) GO TO 550 233 xo1 = xi 234 xi2 = xi * xi 235 Xf1 = a0 + a1 * xi + a2 * xi2 + a3 * xi2 * xi 236 Xprimin = a1 + 2.* a2 * xi + 3.* a3 * xi2 169 237 500 CONTINUE 170 PRINT *,' *** PAS DE SOLUTION **** ',i,xlon2171 STOP 4238 WRITE(6,*) ' Pas de solution ***** ',i,xlon2,iter 239 STOP 6 172 240 550 CONTINUE 173 241 174 xxprim(i) = depi/( FLOAT(iim) * Xprimin) 175 xvrai(i) = depi * (xi - 0.5) + xzoom 176 242 xxprim(i) = depi/ ( FLOAT(iim) * Xprimin ) 243 xvrai(i) = xi + xzoom 177 244 178 245 1500 CONTINUE 179 246 180 247 DO i = 1 , iim 181 xlon (i)= xvrai(i)248 xlon(i) = xvrai(i) 182 249 xprimm(i) = xxprim(i) 183 cc xxlon(i) = xlon(i)*180./pi184 250 ENDDO 185 251 186 cc PRINT *,' XLON avant '187 cc PRINT 68,(xxlon(i),i=1,iim)188 189 252 DO i = 1, iim -1 190 253 IF( xvrai(i+1). LT. xvrai(i) ) THEN 191 PRINT *,' PBS. avec rlonu(',i+1,'plus petit que rlonu(',i,192 , 193 STOP254 WRITE(6,*) ' PBS. avec rlonu(',i+1,') plus petit que rlonu(',i, 255 , ')' 256 STOP 7 194 257 ENDIF 195 258 ENDDO … … 205 268 ENDDO 206 269 207 PRINT *,' LONGITUDES min max ',champmin,champmax208 209 270 IF(champmin .GE. - pi .AND. champmax.LE. pi ) THEN 210 271 GO TO 1600 211 272 ELSE 212 PRINT 18213 PRINT *,'Reorganisation des longitudes pour avoir entre - pi',214 , ' et 273 WRITE(6,18) 274 WRITE(6,*) 'Reorganisation des longitudes pour avoir entre - pi', 275 , ' et pi ' 215 276 c 216 277 IF( xzoom.LE.0.) THEN … … 219 280 IF( xvrai(i).GE. - pi ) GO TO 80 220 281 ENDDO 221 PRINT *, ' PBS. 1'222 STOP 282 WRITE(6,*) ' PBS. 1 ! Xvrai plus petit que - pi ! ' 283 STOP 8 223 284 80 CONTINUE 224 285 is2 = i … … 240 301 IF( xvrai(i).LE. pi ) GO TO 90 241 302 ENDDO 242 PRINT *,' PBS. 2'243 STOP 303 WRITE(6,*) ' PBS. 2 ! Xvrai plus grand que pi ! ' 304 STOP 9 244 305 90 CONTINUE 245 306 is2 = i 246 307 ENDIF 247 cc PRINT *,' IS2 ',is2248 308 idif = iim -is2 249 309 DO ii = 1, is2 … … 270 330 ENDDO 271 331 272 273 332 IF( ik.EQ.1 ) THEN 274 PRINT *, ' XLON aux pts. V-0.25 apres ( en deg. ) ' 275 PRINT 18 276 PRINT 68,xvrai 277 PRINT *,' XPRIM ' 278 PRINT 68, xprimm 333 WRITE(6,*) ' XLON aux pts. V-0.25 apres ( en deg. ) ' 334 WRITE(6,18) 335 WRITE(6,68) xvrai 336 ccc WRITE(6,*) ' XPRIM k ',ik 337 ccc WRITE(6,566) xprimm 338 279 339 DO i = 1,iim + 1 280 340 rlonm025(i) = xlon( i ) 281 341 xprimm025(i) = xprimm(i) 282 342 ENDDO 343 283 344 ELSE IF( ik.EQ.2 ) THEN 284 PRINT 18 285 PRINT *, ' XLON aux pts. V apres ( en deg. ) ' 286 PRINT 68,xvrai 287 PRINT *,' XPRIM ' 288 PRINT 68, xprimm 345 WRITE(6,18) 346 WRITE(6,*) ' XLON aux pts. V apres ( en deg. ) ' 347 WRITE(6,68) xvrai 348 ccc WRITE(6,*) ' XPRIM k ',ik 349 ccc WRITE(6,566) xprimm 350 289 351 DO i = 1,iim + 1 290 352 rlonv(i) = xlon( i ) 291 353 xprimv(i) = xprimm(i) 292 354 ENDDO 293 ELSE IF( ik.EQ.3 ) THEN 294 PRINT 18 295 PRINT *, ' XLON aux pts. U apres ( en deg. ) ' 296 PRINT 68,xvrai 297 PRINT *,' XPRIM ' 298 PRINT 68, xprimm 355 356 ELSE IF( ik.EQ.3) THEN 357 WRITE(6,18) 358 WRITE(6,*) ' XLON aux pts. U apres ( en deg. ) ' 359 WRITE(6,68) xvrai 360 ccc WRITE(6,*) ' XPRIM ik ',ik 361 ccc WRITE(6,566) xprimm 362 299 363 DO i = 1,iim + 1 300 364 rlonu(i) = xlon( i ) 301 365 xprimu(i) = xprimm(i) 302 366 ENDDO 367 303 368 ELSE IF( ik.EQ.4 ) THEN 304 PRINT 18 305 PRINT *, ' XLON aux pts. V+0.25 apres ( en deg. ) ' 306 PRINT 68,xvrai 307 PRINT *,' XPRIM ' 308 PRINT 68, xprimm 369 WRITE(6,18) 370 WRITE(6,*) ' XLON aux pts. V+0.25 apres ( en deg. ) ' 371 WRITE(6,68) xvrai 372 ccc WRITE(6,*) ' XPRIM ik ',ik 373 ccc WRITE(6,566) xprimm 374 309 375 DO i = 1,iim + 1 310 376 rlonp025(i) = xlon( i ) 311 377 xprimp025(i) = xprimm(i) 312 378 ENDDO 379 313 380 ENDIF 314 381 315 382 5000 CONTINUE 316 383 c 317 c ........... fin de la boucle do 5000 ............ 318 319 c 320 DO i = 1, iim + 1 321 dlon1(i) = rlonm025(i) - rlonv(i) 322 dlon2(i) = rlonm025(i) - rlonp025(i) 323 dlon3(i) = rlonm025(i) - rlonu(i) 384 c ........... fin de la boucle do 5000 ............ 385 386 DO i = 1, iim 387 xlon(i) = rlonv(i+1) - rlonv(i) 324 388 ENDDO 325 326 DO i = 1, iim + 1 327 rlonm025(i) = rlonm025(i) + dlon1(i) 389 champmin = 1.e12 390 champmax = -1.e12 391 DO i = 1, iim 392 champmin = MIN( champmin, xlon(i) ) 393 champmax = MAX( champmax, xlon(i) ) 328 394 ENDDO 329 330 DO i = 1, iim + 1 331 rlonv(i) = rlonm025(i) - dlon1(i) 332 rlonp025(i) = rlonm025(i) - dlon2(i) 333 rlonu(i) = rlonm025(i) - dlon3(i) 334 ENDDO 335 336 DO i = 1, iim 337 xprimu (i) = rlonu(i+1) - rlonu(i) 338 xprimv (i) = rlonv(i+1) - rlonv(i) 339 xprimm025(i) = rlonm025(i+1) - rlonm025(i) 340 xprimp025(i) = rlonp025(i+1) - rlonp025(i) 341 ENDDO 342 xprimu (iip1) = xprimu (1) 343 xprimv (iip1) = xprimv (1) 344 xprimm025(iip1) = xprimm025(1) 345 xprimp025(iip1) = xprimp025(1) 346 347 348 18 FORMAT(/) 349 68 FORMAT(1x,7f9.2) 350 351 352 RETURN 353 END 395 champmin = champmin * 180./pi 396 champmax = champmax * 180./pi 397 398 WRITE(6,18) 399 WRITE(6,*) ' Longitudes ' 400 WRITE(6,18) 401 WRITE(6,3) champmin, champmax 402 WRITE(6,*) ' Si cette derniere est trop lache , modifiez les par 403 ,ametres grossism , tau , dzoom pour X et repasser ! ' 404 WRITE(6,18) 405 406 3 Format(1x, ' Au centre du zoom , la longueur de la maille est', 407 , ' d environ ',f8.2 ,' degres ', 408 , ' alors que la maille en dehors de la zone du zoom est d environ 409 , ', f8.2,' degres ' ) 410 18 FORMAT(/) 411 24 FORMAT(2x,' Parametres yzoom,gross,tau ,dzoom pour fyhyp ',4f8.3) 412 68 FORMAT(1x,7f9.2) 413 414 RETURN 415 END -
LMDZ.3.3/trunk/libf/dyn3d/fxyhyper.F
r2 r203 1 SUBROUTINE fxyhyper ( yzoom, grossy, dzoomy,tauy,deltay , 2 , xzoom, grossx, dzoomx,taux, 1 c 2 c $Header 3 c 4 SUBROUTINE fxyhyper ( yzoom, grossy, dzoomy,tauy , 5 , xzoom, grossx, dzoomx,taux , 3 6 , rlatu,yprimu,rlatv,yprimv,rlatu1, yprimu1, rlatu2, yprimu2 , 4 7 , rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025) … … 19 22 c a) le grossissement du zoom : grossy ( en y ) et grossx ( en x ) 20 23 c b) l' extension du zoom : dzoomy ( en y ) et dzoomx ( en x ) 21 c c) la raideur du zoom : tau , ici = 1.24 c c) la raideur de la transition du zoom : taux et tauy 22 25 c 23 c N.B : le produit ( grossissement * extension ) doit etre < 1. 24 c ******* 25 c En plus , il y a un autre parametre , moins important mais quand 26 c meme utile , c'est deltay , deplacement de la zone du zoom en 27 c latitude : Si on deplace de 10. degres vers le nord ( deltay 28 c = 10. dans inigeom ) . 29 c 26 c N.B : Il vaut mieux avoir : grossx * dzoomx < pi ( radians ) 27 c ****** 28 c et grossy * dzoomy < pi/2 ( radians ) 30 29 c 31 30 #include "dimensions.h" … … 40 39 REAL rlonu(iip1),xprimu(iip1),rlonv(iip1),xprimv(iip1), 41 40 , rlonm025(iip1),xprimm025(iip1), rlonp025(iip1),xprimp025(iip1) 42 REAL deltay43 41 44 42 c .... var. locales ..... … … 47 45 c 48 46 49 CALL fyhyp ( yzoom, grossy, dzoomy,tauy , deltay,47 CALL fyhyp ( yzoom, grossy, dzoomy,tauy , 50 48 , rlatu, yprimu,rlatv,yprimv,rlatu2,yprimu2,rlatu1,yprimu1 ) 51 52 49 53 50 CALL fxhyp(xzoom,grossx,dzoomx,taux,rlonm025,xprimm025,rlonv, … … 55 52 56 53 57 DO i = 1, ii m54 DO i = 1, iip1 58 55 IF(rlonp025(i).LT.rlonv(i)) THEN 59 PRINT *,' Attention ! rlonp025 < rlonv',i56 WRITE(6,*) ' Attention ! rlonp025 < rlonv',i 60 57 STOP 61 58 ENDIF 62 59 63 60 IF(rlonv(i).LT.rlonm025(i)) THEN 64 PRINT *,' Attention ! rlonm025 > rlonv',i61 WRITE(6,*) ' Attention ! rlonm025 > rlonv',i 65 62 STOP 66 63 ENDIF 67 64 68 65 IF(rlonp025(i).GT.rlonu(i)) THEN 69 PRINT *,' Attention ! rlonp025 > rlonu',i66 WRITE(6,*) ' Attention ! rlonp025 > rlonu',i 70 67 STOP 71 68 ENDIF 72 69 ENDDO 73 70 74 PRINT *,' *** TEST DE COHERENCE OK POUR FX **** '71 WRITE(6,*) ' *** TEST DE COHERENCE OK POUR FX **** ' 75 72 76 73 c … … 78 75 c 79 76 IF(rlatu1(j).LE.rlatu2(j)) THEN 80 PRINT *,'Attention ! rlatu1 < rlatu2 ',rlatu1(j), rlatu2(j),j77 WRITE(6,*)'Attention ! rlatu1 < rlatu2 ',rlatu1(j), rlatu2(j),j 81 78 STOP 13 82 79 ENDIF 83 80 c 84 81 IF(rlatu2(j).LE.rlatu(j+1)) THEN 85 PRINT *,'Attention ! rlatu2 < rlatup1 ',rlatu2(j),rlatu(j+1),j82 WRITE(6,*)'Attention ! rlatu2 < rlatup1 ',rlatu2(j),rlatu(j+1),j 86 83 STOP 14 87 84 ENDIF 88 85 c 89 86 IF(rlatu(j).LE.rlatu1(j)) THEN 90 PRINT *,' Attention ! rlatu < rlatu1 ',rlatu(j),rlatu1(j),j87 WRITE(6,*)' Attention ! rlatu < rlatu1 ',rlatu(j),rlatu1(j),j 91 88 STOP 15 92 89 ENDIF 93 90 c 94 91 IF(rlatv(j).LE.rlatu2(j)) THEN 95 PRINT *,' Attention ! rlatv < rlatu2 ',rlatv(j),rlatu2(j),j92 WRITE(6,*)' Attention ! rlatv < rlatu2 ',rlatv(j),rlatu2(j),j 96 93 STOP 16 97 94 ENDIF 98 95 c 99 96 IF(rlatv(j).ge.rlatu1(j)) THEN 100 PRINT *,' Attention ! rlatv > rlatu1 ',rlatv(j),rlatu1(j),j97 WRITE(6,*)' Attention ! rlatv > rlatu1 ',rlatv(j),rlatu1(j),j 101 98 STOP 17 102 99 ENDIF 103 100 c 104 101 IF(rlatv(j).ge.rlatu(j)) THEN 105 PRINT *,'Attention ! rlatv > rlatu ',rlatv(j),rlatu(j),j102 WRITE(6,*) ' Attention ! rlatv > rlatu ',rlatv(j),rlatu(j),j 106 103 STOP 18 107 104 ENDIF … … 109 106 ENDDO 110 107 c 111 PRINT *,' *** TEST DE COHERENCE OK POUR FY **** ' 112 108 WRITE(6,*) ' *** TEST DE COHERENCE OK POUR FY **** ' 109 WRITE(6,25) 110 25 FORMAT(//) 113 111 c 114 112 -
LMDZ.3.3/trunk/libf/dyn3d/fyhyp.F
r2 r203 1 SUBROUTINE fyhyp ( yzoomdeg, grossism, dzoom,tau,deltay , 1 c 2 c $Header 3 c 4 SUBROUTINE fyhyp ( yzoomdeg, grossism, dzoom,tau , 2 5 , rrlatu,yyprimu,rrlatv,yyprimv,rlatu2,yprimu2,rlatu1,yprimu1 ) 3 6 7 cc ... Version du 01/04/2001 .... 4 8 5 9 IMPLICIT NONE … … 13 17 c 14 18 c grossism etant le grossissement ( = 2 si 2 fois, = 3 si 3 fois , etc) 15 c dzoom etant la distance totale de la zone du zoom 16 c tau la transition , normalement = 1 . 17 18 c N.B : on doit avoir : grossism * dzoom < 1 . 19 c ************** 20 c 21 c Pour Indoex , on a pris : 22 c ******* 23 c grossism = 2.5 , dzoom = 7/24 en x et y , pour iim = 128 et jjm=64 24 c yzoomdeg = 0. , tau = 1. et delaty = 10. 19 c dzoom etant la distance totale de la zone du zoom ( en radians ) 20 c tau la raideur de la transition de l'interieur a l'exterieur du zoom 21 c 22 c 23 c N.B : Il vaut mieux avoir : grossism * dzoom < pi/2 (radians) ,en lati. 24 c ******************************************************************** 25 25 c 26 26 c … … 29 29 30 30 INTEGER nmax , nmax2 31 PARAMETER ( nmax = 50000, nmax2 = 2*nmax )31 PARAMETER ( nmax = 30000, nmax2 = 2*nmax ) 32 32 c 33 33 c 34 34 c ....... arguments d'entree ....... 35 35 c 36 REAL yzoomdeg, grossism,dzoom,tau , deltay 36 REAL yzoomdeg, grossism,dzoom,tau 37 c ( rentres par run.def ) 37 38 38 39 c ....... arguments de sortie ....... … … 42 43 43 44 c 44 c ..... Champs locaux .....45 c ..... champs locaux ..... 45 46 c 46 47 47 REAl ylat(jjp1), yprim(jjp1) 48 REAL yuv 49 REAL ytild(0:nmax2) 50 REAL fhyp(0:nmax),ffdx(0:nmax),beta,Ytprim(0:nmax2) 51 SAVE Ytprim, ytild,Yf 52 REAL Yf(0:nmax2),yypr(0:nmax2) 53 REAL yvrai(jjp1), yprimm(jjp1),ylatt(jjp1) 54 REAL pi,depi,pis2,epsilon,yzoom 55 REAL yo1,yi,ylon2,fxm,ymoy,yint,Yprimin 56 REAL ypn,deply,y00 48 REAL*8 ylat(jjp1), yprim(jjp1) 49 REAL*8 yuv 50 REAL*8 yt(0:nmax2) 51 REAL*8 fhyp(0:nmax2),beta,Ytprim(0:nmax2),fxm(0:nmax2) 52 SAVE Ytprim, yt,Yf 53 REAL*8 Yf(0:nmax2),yypr(0:nmax2) 54 REAL*8 yvrai(jjp1), yprimm(jjp1),ylatt(jjp1) 55 REAL*8 pi,depi,pis2,epsilon,y0,pisjm 56 REAL*8 yo1,yi,ylon2,ymoy,Yprimin,champmin,champmax 57 REAL*8 yfi,Yf1,ffdy 58 REAL*8 ypn,deply,y00 57 59 SAVE y00, deply 58 60 … … 60 62 INTEGER jpn,jjpn 61 63 SAVE jpn 62 64 REAL*8 a0,a1,a2,a3,yi2,heavyy0,heavyy0m 65 REAL*8 fa(0:nmax2),fb(0:nmax2) 66 REAL y0min,y0max 67 68 REAL*8 heavyside 69 EXTERNAL heavyside 63 70 64 71 pi = 2. * ASIN(1.) 65 72 depi = 2. * pi 66 73 pis2 = pi/2. 67 epsilon = 1.e-6 68 yzoom = yzoomdeg * pi/180. 69 70 74 pisjm = pi/ FLOAT(jjm) 75 epsilon = 1.e-3 76 y0 = yzoomdeg * pi/180. 77 78 IF( dzoom.LT.1.) THEN 79 dzoom = dzoom * pi 80 ELSEIF( dzoom.LT. 12. ) THEN 81 WRITE(6,*) ' Le param. dzoomy pour fyhyp est trop petit ! L aug 82 ,menter et relancer ! ' 83 STOP 1 84 ELSE 85 dzoom = dzoom * pi/180. 86 ENDIF 87 88 WRITE(6,*) ' yzoom ,dzoomy (radians),tau',y0,dzoom,tau 71 89 72 90 DO i = 0, nmax2 73 ytild(i) = FLOAT(i) /nmax2 74 IF( ytild(i).EQ.0.5 ) ytild(i) = ytild(i) + 1.e-6 75 ENDDO 76 77 DO i = 1, nmax 78 fhyp(i) = TANH ( ( ytild(i) - 0.5*(1.- dzoom) ) / 79 , ( tau * ytild(i) * ( 0.5 -ytild(i))) ) 80 ENDDO 81 82 fhyp( 0 ) = - 1. 83 fhyp( nmax ) = 1. 91 yt(i) = - pis2 + FLOAT(i)* pi /nmax2 92 ENDDO 93 94 heavyy0m = heavyside( -y0 ) 95 heavyy0 = heavyside( y0 ) 96 y0min = 2.*y0*heavyy0m - pis2 97 y0max = 2.*y0*heavyy0 + pis2 98 99 DO i = 0, nmax2 100 IF( yt(i).LT.y0 ) THEN 101 fa (i) = tau* (yt(i)-y0+dzoom/2. ) 102 fb(i) = (yt(i)-2.*y0*heavyy0m +pis2) * ( y0 - yt(i) ) 103 ELSEIF ( yt(i).GT.y0 ) THEN 104 fa(i) = tau *(y0-yt(i)+dzoom/2. ) 105 fb(i) = (2.*y0*heavyy0 -yt(i)+pis2) * ( yt(i) - y0 ) 106 ENDIF 107 108 IF( 200.* fb(i) .LT. - fa(i) ) THEN 109 fhyp ( i) = - 1. 110 ELSEIF( 200. * fb(i) .LT. fa(i) ) THEN 111 fhyp ( i) = 1. 112 ELSE 113 fhyp(i) = TANH ( fa(i)/fb(i) ) 114 ENDIF 115 116 IF( yt(i).EQ.y0 ) fhyp(i) = 1. 117 IF(yt(i).EQ. y0min. OR.yt(i).EQ. y0max ) fhyp(i) = -1. 118 119 ENDDO 84 120 85 121 cc .... Calcul de beta .... 86 122 c 87 ffdx( 0 ) = 0. 88 89 DO i = 1, nmax 90 ymoy = 0.5 * ( ytild(i-1) + ytild( i ) ) 91 fxm = TANH ( ( ymoy - 0.5 * ( 1. - dzoom ) ) / 92 , ( tau * ymoy * ( 0.5 -ymoy)) ) 93 ffdx(i) = ffdx(i-1) + fxm * ( ytild(i) - ytild(i-1) ) 94 ENDDO 95 96 beta = ( grossism * ffdx(nmax) - 0.5 ) / ( ffdx(nmax) - 0.5 ) 123 ffdy = 0. 124 125 DO i = 1, nmax2 126 ymoy = 0.5 * ( yt(i-1) + yt( i ) ) 127 IF( ymoy.LT.y0 ) THEN 128 fa(i)= tau * ( ymoy-y0+dzoom/2.) 129 fb(i) = (ymoy-2.*y0*heavyy0m +pis2) * ( y0 - ymoy ) 130 ELSEIF ( ymoy.GT.y0 ) THEN 131 fa(i)= tau * ( y0-ymoy+dzoom/2. ) 132 fb(i) = (2.*y0*heavyy0 -ymoy+pis2) * ( ymoy - y0 ) 133 ENDIF 134 135 IF( 200.* fb(i) .LT. - fa(i) ) THEN 136 fxm ( i) = - 1. 137 ELSEIF( 200. * fb(i) .LT. fa(i) ) THEN 138 fxm ( i) = 1. 139 ELSE 140 fxm(i) = TANH ( fa(i)/fb(i) ) 141 ENDIF 142 IF( ymoy.EQ.y0 ) fxm(i) = 1. 143 IF (ymoy.EQ. y0min. OR.yt(i).EQ. y0max ) fxm(i) = -1. 144 ffdy = ffdy + fxm(i) * ( yt(i) - yt(i-1) ) 145 146 ENDDO 147 148 beta = ( grossism * ffdy - pi ) / ( ffdy - pi ) 149 150 IF( 2.*beta - grossism.LE. 0.) THEN 151 152 WRITE(6,*) ' ** Attention ! La valeur beta calculee dans la rou 153 ,tine fyhyp est mauvaise ! ' 154 WRITE(6,*)'Modifier les valeurs de grossismy ,tauy ou dzoomy', 155 , ' et relancer ! *** ' 156 CALL ABORT 157 158 ENDIF 97 159 c 98 160 c ..... calcul de Ytprim ..... 99 161 c 100 162 101 DO i = 0, nmax 163 DO i = 0, nmax2 102 164 Ytprim(i) = beta + ( grossism - beta ) * fhyp(i) 103 165 ENDDO 104 c105 DO i = 0, nmax106 Ytprim( nmax2 - i ) = Ytprim( i )107 ENDDO108 c109 166 110 167 c ..... Calcul de Yf ........ 111 168 112 Yf(0) = 0. 113 DO i = 1, nmax 114 ymoy = 0.5 * ( ytild(i-1) + ytild( i ) ) 115 fxm = TANH ( ( ymoy - 0.5 * ( 1. - dzoom ) ) / 116 , ( tau * ymoy * ( 0.5 -ymoy)) ) 117 yypr(i) = beta + ( grossism - beta ) * fxm 118 ENDDO 119 120 DO i = 1,nmax 121 yypr(nmax2-i+1) = yypr(i) 122 ENDDO 123 124 DO i=1,nmax2 125 Yf(i) = Yf(i-1) + yypr(i) * ( ytild(i) - ytild(i-1) ) 126 ENDDO 127 c 128 c 169 Yf(0) = - pis2 170 DO i = 1, nmax2 171 yypr(i) = beta + ( grossism - beta ) * fxm(i) 172 ENDDO 173 174 DO i=1,nmax2 175 Yf(i) = Yf(i-1) + yypr(i) * ( yt(i) - yt(i-1) ) 176 ENDDO 129 177 130 178 c **************************************************************** … … 149 197 jlat = jjm 150 198 ENDIF 151 152 c 199 c 200 yo1 = 0. 153 201 DO 1500 j = 1,jlat 154 155 ylon2 = ( FLOAT(j) + yuv -1.) / FLOAT(jjm)156 157 202 yo1 = 0. 158 yi = ylon2 159 160 c 161 DO 500 iter = 1,300 162 203 ylon2 = - pis2 + pisjm * ( FLOAT(j) + yuv -1.) 204 yfi = ylon2 205 c 163 206 DO 250 it = nmax2,0,-1 164 IF( y i.GE.ytild(it)) GO TO 350207 IF( yfi.GE.Yf(it)) GO TO 350 165 208 250 CONTINUE 166 167 209 it = 0 168 yi = ytild(it)169 170 210 350 CONTINUE 171 211 212 yi = yt(it) 172 213 IF(it.EQ.nmax2) THEN 173 214 it = nmax2 -1 174 Yf(it+1) = 1.215 Yf(it+1) = pis2 175 216 ENDIF 176 217 c ................................................................. … … 179 220 c ................................................................. 180 221 181 yint = ( Yf(it+1)-Yf(it) ) / ( ytild(it+1)-ytild(it) ) * 182 + ( yi-ytild(it) ) + Yf(it) 183 Yprimin = ( Ytprim(it+1)-Ytprim(it) )/ ( ytild(it+1)-ytild(it) ) * 184 + ( yi-ytild(it) ) + Ytprim(it) 185 yi = yi - (yint-ylon2)/Yprimin 186 187 IF( ABS(yi-yo1).LE.epsilon) GO TO 550 188 yo1 = yi 189 c 222 CALL coefpoly ( Yf(it),Yf(it+1),Ytprim(it), Ytprim(it+1), 223 , yt(it),yt(it+1) , a0,a1,a2,a3 ) 224 225 Yf1 = Yf(it) 226 Yprimin = a1 + 2.* a2 * yi + 3.*a3 * yi *yi 227 228 DO 500 iter = 1,300 229 yi = yi - ( Yf1 - yfi )/ Yprimin 230 231 IF( ABS(yi-yo1).LE.epsilon) GO TO 550 232 yo1 = yi 233 yi2 = yi * yi 234 Yf1 = a0 + a1 * yi + a2 * yi2 + a3 * yi2 * yi 235 Yprimin = a1 + 2.* a2 * yi + 3.* a3 * yi2 190 236 500 CONTINUE 191 PRINT *,' *** PAS DE SOLUTION ****',j,ylon2,iter192 STOP 4237 WRITE(6,*) ' Pas de solution ***** ',j,ylon2,iter 238 STOP 2 193 239 550 CONTINUE 194 195 yprim(j) = pi /( FLOAT(jjm) * Yprimin) 196 yvrai(j) = pi * (yi - 0.5) + yzoom 240 c 241 Yprimin = a1 + 2.* a2 * yi + 3.* a3 * yi* yi 242 yprim(j) = pi / ( jjm * Yprimin ) 243 yvrai(j) = yi 197 244 198 245 1500 CONTINUE 199 200 cc print *,' LAT avant reorgan '201 cc print 68,(yyvrai(j),j=1,jlat)202 246 203 247 DO j = 1, jlat -1 204 248 IF( yvrai(j+1). LT. yvrai(j) ) THEN 205 PRINT *,' PBS. avec rlat(',j+1,' plus petit que rlat(',j, 206 , ')' 207 STOP 208 ENDIF 209 ENDDO 210 211 PRINT 18 212 PRINT *,'Reorganisation des latitudes pour avoir entre - pi/2 ', 213 , ' et pi/2 ' 214 c 215 249 WRITE(6,*) ' PBS. avec rlat(',j+1,') plus petit que rlat(',j, 250 , ')' 251 STOP 3 252 ENDIF 253 ENDDO 254 255 WRITE(6,18) 256 WRITE(6,*) 'Reorganisation des latitudes pour avoir entre - pi/2' 257 , ,' et pi/2 ' 258 c 216 259 IF( ik.EQ.1 ) THEN 217 ypn = pis2 - deltay * pi/180.260 ypn = pis2 218 261 DO j = jlat,1,-1 219 262 IF( yvrai(j).LE. ypn ) GO TO 1502 … … 239 282 ENDDO 240 283 241 242 243 284 c *********** Fin de la reorganisation ************* 244 285 c 245 286 1600 CONTINUE 246 247 248 287 249 288 DO j = 1, jlat … … 256 295 ENDDO 257 296 258 259 297 IF( ik.EQ.1 ) THEN 260 PRINT 18 261 PRINT *, ' YLAT en U apres ( en deg. ) ' 262 PRINT 68,(yvrai(j),j=1,jlat) 263 PRINT *,' YPRIM ' 264 PRINT 68,( yprim(j),j=1,jlat) 298 WRITE(6,18) 299 WRITE(6,*) ' YLAT en U apres ( en deg. ) ' 300 WRITE(6,68) (yvrai(j),j=1,jlat) 301 cc WRITE(6,*) ' YPRIM ' 302 cc WRITE(6,445) ( yprim(j),j=1,jlat) 303 265 304 DO j = 1, jlat 266 305 rrlatu(j) = ylat( j ) 267 306 yyprimu(j) = yprim( j ) 268 307 ENDDO 269 c 308 270 309 ELSE IF ( ik.EQ. 2 ) THEN 271 PRINT 18 272 PRINT *, ' YLAT en V apres ( en deg. ) ' 273 PRINT 68,(yvrai(j),j=1,jlat) 274 PRINT *,' YPRIM ' 275 PRINT 68,( yprim(j),j=1,jlat) 310 WRITE(6,18) 311 WRITE(6,*) ' YLAT en V apres ( en deg. ) ' 312 WRITE(6,68) (yvrai(j),j=1,jlat) 313 cc WRITE(6,*)' YPRIM ' 314 cc WRITE(6,445) ( yprim(j),j=1,jlat) 315 276 316 DO j = 1, jlat 277 317 rrlatv(j) = ylat( j ) 278 318 yyprimv(j) = yprim( j ) 279 319 ENDDO 280 c 320 281 321 ELSE IF ( ik.EQ. 3 ) THEN 282 PRINT 18 283 PRINT *, ' YLAT en U + 0.75 apres ( en deg. ) ' 284 PRINT 68,(yvrai(j),j=1,jlat) 285 PRINT *,' YPRIM ' 286 PRINT 68,( yprim(j),j=1,jlat) 322 WRITE(6,18) 323 WRITE(6,*) ' YLAT en U + 0.75 apres ( en deg. ) ' 324 WRITE(6,68) (yvrai(j),j=1,jlat) 325 cc WRITE(6,*) ' YPRIM ' 326 cc WRITE(6,445) ( yprim(j),j=1,jlat) 327 287 328 DO j = 1, jlat 288 329 rlatu2(j) = ylat( j ) … … 291 332 292 333 ELSE IF ( ik.EQ. 4 ) THEN 293 PRINT 18 294 PRINT *, ' YLAT en U + 0.25 apres ( en deg. ) ' 295 PRINT 68,(yvrai(j),j=1,jlat) 296 PRINT *,' YPRIM ' 297 PRINT 68,( yprim(j),j=1,jlat) 334 WRITE(6,18) 335 WRITE(6,*) ' YLAT en U + 0.25 apres ( en deg. ) ' 336 WRITE(6,68)(yvrai(j),j=1,jlat) 337 cc WRITE(6,*) ' YPRIM ' 338 cc WRITE(6,68) ( yprim(j),j=1,jlat) 339 298 340 DO j = 1, jlat 299 341 rlatu1(j) = ylat( j ) 300 342 yprimu1(j) = yprim( j ) 301 343 ENDDO 344 302 345 ENDIF 303 346 … … 306 349 c ..... fin de la boucle do 5000 ..... 307 350 351 DO j = 1, jjm 352 ylat(j) = rrlatu(j) - rrlatu(j+1) 353 ENDDO 354 champmin = 1.e12 355 champmax = -1.e12 356 DO j = 1, jjm 357 champmin = MIN( champmin, ylat(j) ) 358 champmax = MAX( champmax, ylat(j) ) 359 ENDDO 360 champmin = champmin * 180./pi 361 champmax = champmax * 180./pi 362 363 WRITE(6,18) 364 WRITE(6,*) ' Latitudes ' 365 WRITE(6,18) 366 WRITE(6,3) champmin, champmax 367 WRITE(6,*) ' Si cette derniere est trop lache , modifiez les par 368 ,ametres grossism , tau , dzoom pour Y et repasser ! ' 369 WRITE(6,18) 370 c 371 3 Format(1x, ' Au centre du zoom , la longueur de la maille est', 372 , ' d environ ',f8.2 ,' degres ', 373 , ' alors que la maille en dehors de la zone du zoom est d environ 374 , ', f8.2,' degres ' ) 308 375 18 FORMAT(/) 309 376 68 FORMAT(1x,7f9.2) -
LMDZ.3.3/trunk/libf/dyn3d/inigeom.F
r2 r203 1 c 2 c $Header 3 c 1 4 SUBROUTINE inigeom 2 5 c 3 6 c Auteur : P. Le Van 4 c ..................... 5 c 6 c ............ Version du 20/12/98 ........................ 7 c 8 c ............ Version du 01/04/2001 ........................ 7 9 c 8 10 c Calcul des elongations cuij1,.cuij4 , cvij1,..cvij4 aux memes en- 9 11 c endroits que les aires aireij1,..aireij4 . 12 10 13 c Choix entre f(y) a derivee sinusoid. ou a derivee tangente hyperbol. 11 14 c 12 15 c 13 16 IMPLICIT NONE 14 c15 REAL deltay, tauy, taux16 PARAMETER ( deltay = 0. , tauy = 1. , taux = 1. )17 c18 c deltay est ( en degres ) le deplacement eventuel en Y du zoom19 cc taux et tauy sont les raideurs du zoom20 c21 17 c 22 18 #include "dimensions.h" … … 50 46 SAVE rlonm025,xprimm025,rlonp025,xprimp025 51 47 52 53 54 c----------------------------------------------------------------------55 48 REAL SSUM 56 49 EXTERNAL SSUM 57 c58 50 c 59 51 c … … 62 54 c - calcul des coeff. ( cu, cv , 1./cu**2, 1./cv**2 ) - 63 55 c - - 64 c ------------------------------------------------------------------65 56 c ------------------------------------------------------------------ 66 57 c … … 168 159 c 169 160 c 170 PRINT 3161 WRITE(6,3) 171 162 3 FORMAT( // 10x,' .... INIGEOM date du 01/06/98 ..... ', 172 163 * //5x,' Calcul des elongations cu et cv comme sommes des 4 ' / … … 191 182 ENDIF 192 183 193 PRINT *,' gamdi_gd ',gamdi_gdiv,gamdi_grot,gamdi_h,coefdis,184 WRITE(6,*) ' gamdi_gd ',gamdi_gdiv,gamdi_grot,gamdi_h,coefdis, 194 185 * nitergdiv,nitergrot,niterh 195 186 c 196 187 pi = 2.* ASIN(1.) 197 188 c 198 PRINT 990189 WRITE(6,990) 199 190 200 191 c ---------------------------------------------------------------- … … 205 196 IF( ysinus ) THEN 206 197 c 207 PRINT *,' *** Inigeom , Y = Sinus ( Latitude ) *** '198 WRITE(6,*) ' *** Inigeom , Y = Sinus ( Latitude ) *** ' 208 199 c 209 200 c .... utilisation de f(x,y ) avec y = sinus de la latitude ..... … … 215 206 ELSE 216 207 c 217 PRINT *,'*** Inigeom , Y = Latitude , der. sinusoid . ***'208 WRITE(6,*) '*** Inigeom , Y = Latitude , der. sinusoid . ***' 218 209 219 210 c .... utilisation de f(x,y) a tangente sinusoidale , y etant la latit. ... … … 267 258 ELSE 268 259 c 269 c .... utilisation de fxyhyper , f(x,y) a derivee tangente hyperbol.260 c .... Utilisation de fxyhyper , f(x,y) a derivee tangente hyperbol. 270 261 c ..................................................................... 271 262 272 PRINT *,'*** Inigeom , Y = Latitude , der.tg. hyperbolique ***'273 274 CALL fxyhyper( clat, grossismy, dzoomy, tauy , deltay ,275 , clon, grossismx, dzoomx, taux,263 WRITE(6,*)'*** Inigeom , Y = Latitude , der.tg. hyperbolique ***' 264 265 CALL fxyhyper( clat, grossismy, dzoomy, tauy , 266 , clon, grossismx, dzoomx, taux , 276 267 , rlatu,yprimu,rlatv, yprimv,rlatu1, yprimu1,rlatu2,yprimu2 , 277 268 , rlonu,xprimu,rlonv,xprimv,rlonm025,xprimm025,rlonp025,xprimp025 ) … … 387 378 c 388 379 IF ( j. eq. jjp1 ) THEN 389 390 380 yprp = yprimu2(j-1) 391 381 rlatp = rlatu2 (j-1) 392 cc yprp = fyprim( FLOAT(j) - 0.25 )393 cc rlatp = fy ( FLOAT(j) - 0.25 )382 ccc yprp = fyprim( FLOAT(j) - 0.25 ) 383 ccc rlatp = fy ( FLOAT(j) - 0.25 ) 394 384 c 395 385 coslatp = COS( rlatp ) … … 425 415 rlatm = rlatu1 ( j ) 426 416 yprm = yprimu1( j ) 427 c rlatp = fy ( FLOAT(j) - 0.25 )428 c yprp = fyprim( FLOAT(j) - 0.25 )429 c rlatm = fy ( FLOAT(j) + 0.25 )430 c yprm = fyprim( FLOAT(j) + 0.25 )417 cc rlatp = fy ( FLOAT(j) - 0.25 ) 418 cc yprp = fyprim( FLOAT(j) - 0.25 ) 419 cc rlatm = fy ( FLOAT(j) + 0.25 ) 420 cc yprm = fyprim( FLOAT(j) + 0.25 ) 431 421 432 422 coslatm = COS( rlatm ) … … 490 480 36 CONTINUE 491 481 c 492 c .... Modif P. Le Van ( 4/07/96 ) .....493 482 c 494 483 aire (iip1,j) = aire (1,j) … … 620 609 aiuscv2gam(iip1,j) = aiuscv2gam(1,j) 621 610 ENDDO 622 c 611 623 612 c 624 613 c calcul des aires aux poles : … … 666 655 c----------------------------------------------------------------------- 667 656 c 668 PRINT *,' INIGEOM RLONV ' 657 WRITE(6,*) ' *** Coordonnees de la grille *** ' 658 WRITE(6,995) 659 c 660 WRITE(6,*) ' LONGITUDES aux pts. V ( degres ) ' 661 WRITE(6,995) 669 662 DO i=1,iip1 670 663 rlonvv(i) = rlonv(i)*180./pi 671 664 ENDDO 672 PRINT 400,rlonvv 673 c 674 PRINT *,' RLATV ' 665 WRITE(6,400) rlonvv 666 c 667 WRITE(6,995) 668 WRITE(6,*) ' LATITUDES aux pts. V ( degres ) ' 669 WRITE(6,995) 675 670 DO i=1,jjm 676 671 rlatuu(i)=rlatv(i)*180./pi 677 672 ENDDO 678 PRINT 400,(rlatuu(i),i=1,jjm)673 WRITE(6,400) (rlatuu(i),i=1,jjm) 679 674 c 680 675 DO i=1,iip1 681 676 rlonvv(i)=rlonu(i)*180./pi 682 677 ENDDO 683 PRINT *,' RLONU ' 684 PRINT 400,rlonvv 685 c 686 PRINT *,' RLATU ' 678 WRITE(6,995) 679 WRITE(6,*) ' LONGITUDES aux pts. U ( degres ) ' 680 WRITE(6,995) 681 WRITE(6,400) rlonvv 682 WRITE(6,995) 683 684 WRITE(6,*) ' LATITUDES aux pts. U ( degres ) ' 685 WRITE(6,995) 687 686 DO i=1,jjp1 688 687 rlatuu(i)=rlatu(i)*180./pi 689 688 ENDDO 690 PRINT 400,(rlatuu(i),i=1,jjp1) 691 c 689 WRITE(6,400) (rlatuu(i),i=1,jjp1) 690 WRITE(6,995) 691 c 692 444 format(f10.3,f6.0) 692 693 400 FORMAT(1x,8f8.2) 693 694 990 FORMAT(//) 695 995 FORMAT(/) 694 696 c 695 697 RETURN -
LMDZ.3.3/trunk/libf/dyn3d/serre.h
r2 r203 1 c 2 c $Header 3 c 1 4 c..include serre.h 2 5 c 3 6 REAL clon,clat,transx,transy,alphax,alphay,pxo,pyo, 4 , grossismx, grossismy, dzoomx, dzoomy 7 , grossismx, grossismy, dzoomx, dzoomy,taux,tauy 5 8 COMMON/serre/clon,clat,transx,transy,alphax,alphay,pxo,pyo , 6 , grossismx, grossismy, dzoomx, dzoomy 9 , grossismx, grossismy, dzoomx, dzoomy,taux,tauy
Note: See TracChangeset
for help on using the changeset viewer.
