Changeset 298 for LMDZ.3.3/branches/rel-1-0-patch/libf/dyn3d/fxhyp.F
- Timestamp:
- Nov 29, 2001, 1:15:59 PM (23 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ.3.3/branches/rel-1-0-patch/libf/dyn3d/fxhyp.F
r253 r298 20 20 INTEGER nmax, nmax2 21 21 PARAMETER ( nmax = 30000, nmax2 = 2*nmax ) 22 c 23 LOGICAL scal180 24 PARAMETER ( scal180 = .TRUE. ) 25 26 c scal180 = .TRUE. si on veut avoir le premier point scalaire pour 27 c une grille reguliere ( grossism = 1.,tau=0.,clon=0. ) a -180. degres. 28 c sinon scal180 = .FALSE. 22 29 23 30 #include "dimensions.h" … … 42 49 REAL*8 pi,depi,epsilon,xzoom,fa,fb 43 50 REAL*8 Xf1, Xfi , a0,a1,a2,a3,xi2 44 INTEGER i,it,ik,iter,ii,idif 51 INTEGER i,it,ik,iter,ii,idif,ii1,ii2 45 52 REAL*8 xi,xo1,xmoy,xlon2,fxm,Xprimin 46 REAL*8 champmin,champmax 53 REAL*8 champmin,champmax,decalx 47 54 INTEGER is2 48 55 SAVE is2 … … 55 62 epsilon = 1.e-3 56 63 xzoom = xzoomdeg * pi/180. 57 64 c 65 decalx = .75 66 IF( grossism.EQ.1..AND.scal180 ) THEN 67 decalx = 1. 68 ENDIF 69 70 WRITE(6,*) 'FXHYP scal180,decalx', scal180,decalx 71 c 58 72 IF( dzoom.LT.1.) THEN 59 73 dzoom = dzoom * depi … … 73 87 ENDDO 74 88 75 DO i = 0, nmax289 DO i = nmax, nmax2 76 90 77 91 fa = tau* ( dzoom/2. - xtild(i) ) 78 92 fb = xtild(i) * ( pi - xtild(i) ) 79 93 80 81 if ( xtild(i) .ne. 0. .and. xtild(i) .ne. pi) then82 94 IF( 200.* fb .LT. - fa ) THEN 83 95 fhyp ( i) = - 1. … … 85 97 fhyp ( i) = 1. 86 98 ELSE 87 fhyp(i) = TANH ( fa/fb ) 99 IF( ABS(fa).LT.1.e-13.AND.ABS(fb).LT.1.e-13) THEN 100 IF( 200.*fb + fa.LT.1.e-10 ) THEN 101 fhyp ( i ) = - 1. 102 ELSEIF( 200.*fb - fa.LT.1.e-10 ) THEN 103 fhyp ( i ) = 1. 104 ENDIF 105 ELSE 106 fhyp ( i ) = TANH ( fa/fb ) 107 ENDIF 88 108 ENDIF 89 endif 90 91 IF ( xtild(i).EQ. 0. ) fhyp(i) = 1. 92 IF ( xtild(i).EQ. pi ) fhyp(i) = -1. 109 110 IF ( xtild(i).EQ. 0. ) fhyp(i) = 1. 111 IF ( xtild(i).EQ. pi ) fhyp(i) = -1. 93 112 94 113 ENDDO … … 100 119 101 120 DO i = nmax +1,nmax2 121 122 xmoy = 0.5 * ( xtild(i-1) + xtild( i ) ) 123 fa = tau* ( dzoom/2. - xmoy ) 124 fb = xmoy * ( pi - xmoy ) 125 126 IF( 200.* fb .LT. - fa ) THEN 127 fxm = - 1. 128 ELSEIF( 200. * fb .LT. fa ) THEN 129 fxm = 1. 130 ELSE 131 IF( ABS(fa).LT.1.e-13.AND.ABS(fb).LT.1.e-13) THEN 132 IF( 200.*fb + fa.LT.1.e-10 ) THEN 133 fxm = - 1. 134 ELSEIF( 200.*fb - fa.LT.1.e-10 ) THEN 135 fxm = 1. 136 ENDIF 137 ELSE 138 fxm = TANH ( fa/fb ) 139 ENDIF 140 ENDIF 141 142 IF ( xmoy.EQ. 0. ) fxm = 1. 143 IF ( xmoy.EQ. pi ) fxm = -1. 144 145 ffdx = ffdx + fxm * ( xtild(i) - xtild(i-1) ) 146 147 ENDDO 148 149 beta = ( grossism * ffdx - pi ) / ( ffdx - pi ) 150 151 IF( 2.*beta - grossism.LE. 0.) THEN 152 WRITE(6,*) ' ** Attention ! La valeur beta calculee dans la rou 153 ,tine fxhyp est mauvaise ! ' 154 WRITE(6,*)'Modifier les valeurs de grossismx ,tau ou dzoomx ', 155 , ' et relancer ! *** ' 156 CALL ABORT 157 ENDIF 158 c 159 c ..... calcul de Xprimt ..... 160 c 161 162 DO i = nmax, nmax2 163 Xprimt(i) = beta + ( grossism - beta ) * fhyp(i) 164 ENDDO 165 c 166 DO i = nmax+1, nmax2 167 Xprimt( nmax2 - i ) = Xprimt( i ) 168 ENDDO 169 c 170 171 c ..... Calcul de Xf ........ 172 173 Xf(0) = - pi 174 175 DO i = nmax +1, nmax2 102 176 103 177 xmoy = 0.5 * ( xtild(i-1) + xtild( i ) ) … … 113 187 ENDIF 114 188 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 189 IF ( xmoy.EQ. 0. ) fxm = 1. 162 190 IF ( xmoy.EQ. pi ) fxm = -1. … … 185 213 186 214 IF( ik.EQ.1 ) THEN 187 xuv = -0.25215 xuv = -0.25 188 216 ELSE IF ( ik.EQ.2 ) THEN 189 217 xuv = 0. 190 218 ELSE IF ( ik.EQ.3 ) THEN 191 xuv = 0.5 219 xuv = 0.50 192 220 ELSE IF ( ik.EQ.4 ) THEN 193 221 xuv = 0.25 … … 196 224 xo1 = 0. 197 225 198 DO 1500 i = 1, iim 199 200 xlon2 = - pi + (FLOAT(i) + xuv - 0.75) * depi / FLOAT(iim) 226 ii1=1 227 ii2=iim 228 IF(ik.EQ.1.and.grossism.EQ.1.) THEN 229 ii1 = 2 230 ii2 = iim+1 231 ENDIF 232 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 294 255 295 DO i = 1, iim -1 256 296 IF( xvrai(i+1). LT. xvrai(i) ) THEN … … 271 311 ENDDO 272 312 273 IF(champmin .GE. - pi .AND. champmax.LE. pi) THEN313 IF(champmin .GE.-pi-0.10.and.champmax.LE.pi+0.10 ) THEN 274 314 GO TO 1600 275 315 ELSE … … 301 341 IF( ik.EQ.1 ) THEN 302 342 DO i = iim,1,-1 303 IF( xvrai(i).LE. pi )GO TO 90343 IF( xvrai(i).LE. pi ) GO TO 90 304 344 ENDDO 305 345 WRITE(6,*) ' PBS. 2 ! Xvrai plus grand que pi ! ' … … 336 376 c WRITE(6,18) 337 377 c WRITE(6,68) xvrai 338 c ccWRITE(6,*) ' XPRIM k ',ik339 c ccWRITE(6,566) xprimm340 341 DO i = 1,iim + 378 c WRITE(6,*) ' XPRIM k ',ik 379 c WRITE(6,566) xprimm 380 381 DO i = 1,iim +1 342 382 rlonm025(i) = xlon( i ) 343 383 xprimm025(i) = xprimm(i) 344 384 ENDDO 345 385 346 386 ELSE IF( ik.EQ.2 ) THEN 347 387 c WRITE(6,18) 348 388 c WRITE(6,*) ' XLON aux pts. V apres ( en deg. ) ' 349 389 c WRITE(6,68) xvrai 350 c ccWRITE(6,*) ' XPRIM k ',ik351 c ccWRITE(6,566) xprimm390 c WRITE(6,*) ' XPRIM k ',ik 391 c WRITE(6,566) xprimm 352 392 353 393 DO i = 1,iim + 1 … … 360 400 c WRITE(6,*) ' XLON aux pts. U apres ( en deg. ) ' 361 401 c WRITE(6,68) xvrai 362 c ccWRITE(6,*) ' XPRIM ik ',ik363 c ccWRITE(6,566) xprimm402 c WRITE(6,*) ' XPRIM ik ',ik 403 c WRITE(6,566) xprimm 364 404 365 405 DO i = 1,iim + 1 … … 372 412 c WRITE(6,*) ' XLON aux pts. V+0.25 apres ( en deg. ) ' 373 413 c WRITE(6,68) xvrai 374 c ccWRITE(6,*) ' XPRIM ik ',ik375 c ccWRITE(6,566) xprimm414 c WRITE(6,*) ' XPRIM ik ',ik 415 c WRITE(6,566) xprimm 376 416 377 417 DO i = 1,iim + 1 … … 403 443 24 FORMAT(2x,'Parametres xzoom,gross,tau ,dzoom pour fxhyp ',4f8.3) 404 444 68 FORMAT(1x,7f9.2) 445 566 FORMAT(1x,7f9.4) 405 446 406 447 RETURN
Note: See TracChangeset
for help on using the changeset viewer.
