- Timestamp:
- Feb 6, 2002, 3:53:54 PM (23 years ago)
- File:
-
- 1 edited
-
LMDZ.3.3/branches/rel-LF/libf/dyn3d/fxhyp.F (modified) (16 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ.3.3/branches/rel-LF/libf/dyn3d/fxhyp.F
r242 r330 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, … … 20 23 INTEGER nmax, nmax2 21 24 PARAMETER ( nmax = 30000, nmax2 = 2*nmax ) 25 c 26 LOGICAL scal180 27 PARAMETER ( scal180 = .TRUE. ) 28 29 c scal180 = .TRUE. si on veut avoir le premier point scalaire pour 30 c une grille reguliere ( grossism = 1.,tau=0.,clon=0. ) a -180. degres. 31 c sinon scal180 = .FALSE. 22 32 23 33 #include "dimensions.h" … … 42 52 REAL*8 pi,depi,epsilon,xzoom,fa,fb 43 53 REAL*8 Xf1, Xfi , a0,a1,a2,a3,xi2 44 INTEGER i,it,ik,iter,ii,idif 54 INTEGER i,it,ik,iter,ii,idif,ii1,ii2 45 55 REAL*8 xi,xo1,xmoy,xlon2,fxm,Xprimin 46 REAL*8 champmin,champmax 56 REAL*8 champmin,champmax,decalx 47 57 INTEGER is2 48 58 SAVE is2 … … 55 65 epsilon = 1.e-3 56 66 xzoom = xzoomdeg * pi/180. 57 67 c 68 decalx = .75 69 IF( grossism.EQ.1..AND.scal180 ) THEN 70 decalx = 1. 71 ENDIF 72 73 WRITE(6,*) 'FXHYP scal180,decalx', scal180,decalx 74 c 58 75 IF( dzoom.LT.1.) THEN 59 76 dzoom = dzoom * depi … … 73 90 ENDDO 74 91 75 DO i = 0, nmax292 DO i = nmax, nmax2 76 93 77 94 fa = tau* ( dzoom/2. - xtild(i) ) 78 95 fb = xtild(i) * ( pi - xtild(i) ) 79 96 80 81 if ( xtild(i) .ne. 0. .and. xtild(i) .ne. pi) then82 97 IF( 200.* fb .LT. - fa ) THEN 83 98 fhyp ( i) = - 1. … … 85 100 fhyp ( i) = 1. 86 101 ELSE 87 fhyp(i) = TANH ( fa/fb ) 102 IF( ABS(fa).LT.1.e-13.AND.ABS(fb).LT.1.e-13) THEN 103 IF( 200.*fb + fa.LT.1.e-10 ) THEN 104 fhyp ( i ) = - 1. 105 ELSEIF( 200.*fb - fa.LT.1.e-10 ) THEN 106 fhyp ( i ) = 1. 107 ENDIF 108 ELSE 109 fhyp ( i ) = TANH ( fa/fb ) 110 ENDIF 88 111 ENDIF 89 endif 90 91 IF ( xtild(i).EQ. 0. ) fhyp(i) = 1. 92 IF ( xtild(i).EQ. pi ) fhyp(i) = -1. 112 IF ( xtild(i).EQ. 0. ) fhyp(i) = 1. 113 IF ( xtild(i).EQ. pi ) fhyp(i) = -1. 93 114 94 115 ENDDO 95 116 96 117 cc .... Calcul de beta .... 97 c ............................98 118 99 119 ffdx = 0. 100 120 101 121 DO i = nmax +1,nmax2 122 123 xmoy = 0.5 * ( xtild(i-1) + xtild( i ) ) 124 fa = tau* ( dzoom/2. - xmoy ) 125 fb = xmoy * ( pi - xmoy ) 126 127 IF( 200.* fb .LT. - fa ) THEN 128 fxm = - 1. 129 ELSEIF( 200. * fb .LT. fa ) THEN 130 fxm = 1. 131 ELSE 132 IF( ABS(fa).LT.1.e-13.AND.ABS(fb).LT.1.e-13) THEN 133 IF( 200.*fb + fa.LT.1.e-10 ) THEN 134 fxm = - 1. 135 ELSEIF( 200.*fb - fa.LT.1.e-10 ) THEN 136 fxm = 1. 137 ENDIF 138 ELSE 139 fxm = TANH ( fa/fb ) 140 ENDIF 141 ENDIF 142 143 IF ( xmoy.EQ. 0. ) fxm = 1. 144 IF ( xmoy.EQ. pi ) fxm = -1. 145 146 ffdx = ffdx + fxm * ( xtild(i) - xtild(i-1) ) 147 148 ENDDO 149 150 beta = ( grossism * ffdx - pi ) / ( ffdx - pi ) 151 152 IF( 2.*beta - grossism.LE. 0.) THEN 153 WRITE(6,*) ' ** Attention ! La valeur beta calculee dans la rou 154 ,tine fxhyp est mauvaise ! ' 155 WRITE(6,*)'Modifier les valeurs de grossismx ,tau ou dzoomx ', 156 , ' et relancer ! *** ' 157 CALL ABORT 158 ENDIF 159 c 160 c ..... calcul de Xprimt ..... 161 c 162 163 DO i = nmax, nmax2 164 Xprimt(i) = beta + ( grossism - beta ) * fhyp(i) 165 ENDDO 166 c 167 DO i = nmax+1, nmax2 168 Xprimt( nmax2 - i ) = Xprimt( i ) 169 ENDDO 170 c 171 172 c ..... Calcul de Xf ........ 173 174 Xf(0) = - pi 175 176 DO i = nmax +1, nmax2 102 177 103 178 xmoy = 0.5 * ( xtild(i-1) + xtild( i ) ) … … 113 188 ENDIF 114 189 115 IF ( xmoy.EQ. 0. ) fhyp(i) = 1.116 IF ( xmoy.EQ. pi ) fhyp(i) = -1.117 ffdx = ffdx + fxm * ( xtild(i) - xtild(i-1) )118 119 ENDDO120 121 beta = ( grossism * ffdx - pi ) / ( ffdx - pi )122 123 IF( 2.*beta - grossism.LE. 0.) THEN124 WRITE(6,*) ' ** Attention ! La valeur beta calculee dans la rou125 ,tine fxhyp est mauvaise ! '126 WRITE(6,*)'Modifier les valeurs de grossismx ,tau ou dzoomx ',127 , ' et relancer ! *** '128 CALL ABORT129 ENDIF130 c131 c ..... calcul de Xprimt .....132 c133 134 DO i = nmax, nmax2135 Xprimt(i) = beta + ( grossism - beta ) * fhyp(i)136 ENDDO137 c138 DO i = nmax+1, nmax2139 Xprimt( nmax2 - i ) = Xprimt( i )140 ENDDO141 c142 143 c ..... Calcul de Xf ........144 145 Xf(0) = - pi146 147 DO i = nmax +1, nmax2148 149 xmoy = 0.5 * ( xtild(i-1) + xtild( i ) )150 fa = tau* ( dzoom/2. - xmoy )151 fb = xmoy * ( pi - xmoy )152 153 IF( 200.* fb .LT. - fa ) THEN154 fxm = - 1.155 ELSEIF( 200. * fb .LT. fa ) THEN156 fxm = 1.157 ELSE158 fxm = TANH ( fa/fb )159 ENDIF160 161 190 IF ( xmoy.EQ. 0. ) fxm = 1. 162 191 IF ( xmoy.EQ. pi ) fxm = -1. … … 185 214 186 215 IF( ik.EQ.1 ) THEN 187 xuv = -0.25216 xuv = -0.25 188 217 ELSE IF ( ik.EQ.2 ) THEN 189 218 xuv = 0. 190 219 ELSE IF ( ik.EQ.3 ) THEN 191 xuv = 0.5 220 xuv = 0.50 192 221 ELSE IF ( ik.EQ.4 ) THEN 193 222 xuv = 0.25 … … 196 225 xo1 = 0. 197 226 198 DO 1500 i = 1, iim 199 200 xlon2 = - pi + (FLOAT(i) + xuv - 0.75) * depi / FLOAT(iim) 227 ii1=1 228 ii2=iim 229 IF(ik.EQ.1.and.grossism.EQ.1.) THEN 230 ii1 = 2 231 ii2 = iim+1 232 ENDIF 233 DO 1500 i = ii1, ii2 234 235 xlon2 = - pi + (FLOAT(i) + xuv - decalx) * depi / FLOAT(iim) 201 236 202 237 Xfi = xlon2 … … 248 283 1500 CONTINUE 249 284 285 286 IF(ik.EQ.1.and.grossism.EQ.1.) THEN 287 xvrai(1) = xvrai(iip1)-depi 288 xxprim(1) = xxprim(iip1) 289 ENDIF 250 290 DO i = 1 , iim 251 291 xlon(i) = xvrai(i) 252 292 xprimm(i) = xxprim(i) 253 293 ENDDO 254 255 294 DO i = 1, iim -1 256 295 IF( xvrai(i+1). LT. xvrai(i) ) THEN … … 271 310 ENDDO 272 311 273 IF(champmin .GE. - pi .AND. champmax.LE. pi) THEN312 IF(champmin .GE.-pi-0.10.and.champmax.LE.pi+0.10 ) THEN 274 313 GO TO 1600 275 314 ELSE … … 301 340 IF( ik.EQ.1 ) THEN 302 341 DO i = iim,1,-1 303 IF( xvrai(i).LE. pi )GO TO 90342 IF( xvrai(i).LE. pi ) GO TO 90 304 343 ENDDO 305 344 WRITE(6,*) ' PBS. 2 ! Xvrai plus grand que pi ! ' … … 336 375 c WRITE(6,18) 337 376 c WRITE(6,68) xvrai 338 c ccWRITE(6,*) ' XPRIM k ',ik339 c ccWRITE(6,566) xprimm340 341 DO i = 1,iim + 377 c WRITE(6,*) ' XPRIM k ',ik 378 c WRITE(6,566) xprimm 379 380 DO i = 1,iim +1 342 381 rlonm025(i) = xlon( i ) 343 382 xprimm025(i) = xprimm(i) 344 383 ENDDO 345 346 384 ELSE IF( ik.EQ.2 ) THEN 347 385 c WRITE(6,18) 348 386 c WRITE(6,*) ' XLON aux pts. V apres ( en deg. ) ' 349 387 c WRITE(6,68) xvrai 350 c ccWRITE(6,*) ' XPRIM k ',ik351 c ccWRITE(6,566) xprimm388 c WRITE(6,*) ' XPRIM k ',ik 389 c WRITE(6,566) xprimm 352 390 353 391 DO i = 1,iim + 1 … … 360 398 c WRITE(6,*) ' XLON aux pts. U apres ( en deg. ) ' 361 399 c WRITE(6,68) xvrai 362 c ccWRITE(6,*) ' XPRIM ik ',ik363 c ccWRITE(6,566) xprimm400 c WRITE(6,*) ' XPRIM ik ',ik 401 c WRITE(6,566) xprimm 364 402 365 403 DO i = 1,iim + 1 … … 372 410 c WRITE(6,*) ' XLON aux pts. V+0.25 apres ( en deg. ) ' 373 411 c WRITE(6,68) xvrai 374 c ccWRITE(6,*) ' XPRIM ik ',ik375 c ccWRITE(6,566) xprimm412 c WRITE(6,*) ' XPRIM ik ',ik 413 c WRITE(6,566) xprimm 376 414 377 415 DO i = 1,iim + 1 … … 403 441 24 FORMAT(2x,'Parametres xzoom,gross,tau ,dzoom pour fxhyp ',4f8.3) 404 442 68 FORMAT(1x,7f9.2) 443 566 FORMAT(1x,7f9.4) 405 444 406 445 RETURN
Note: See TracChangeset
for help on using the changeset viewer.
