Ignore:
Timestamp:
Nov 29, 2001, 1:09:18 PM (23 years ago)
Author:
lmdz
Message:

Probleme de decalage de grille par rapport a 180 degres regle (Levan)
LF

File:
1 edited

Legend:

Unmodified
Added
Removed
  • LMDZ.3.3/trunk/libf/dyn3d/fxhyp.F

    r251 r297  
     1c
     2c $Header$
     3c
    14       SUBROUTINE fxhyp ( xzoomdeg,grossism,dzoom,tau ,
    25     , rlonm025,xprimm025,rlonv,xprimv,rlonu,xprimu,rlonp025,xprimp025,
     
    2023       INTEGER nmax, nmax2
    2124       PARAMETER (  nmax = 30000, nmax2 = 2*nmax )
     25c
     26       LOGICAL scal180
     27       PARAMETER ( scal180 = .TRUE. )
     28
     29c      scal180 = .TRUE.  si on veut avoir le premier point scalaire pour   
     30c      une grille reguliere ( grossism = 1.,tau=0.,clon=0. ) a -180. degres.
     31c      sinon scal180 = .FALSE.
    2232
    2333#include "dimensions.h"
     
    4252       REAL*8 pi,depi,epsilon,xzoom,fa,fb
    4353       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
    4555       REAL*8 xi,xo1,xmoy,xlon2,fxm,Xprimin
    46        REAL*8 champmin,champmax
     56       REAL*8 champmin,champmax,decalx
    4757       INTEGER is2
    4858       SAVE is2
     
    5565       epsilon  = 1.e-3
    5666       xzoom    = xzoomdeg * pi/180.
    57 
     67c
     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
     74c
    5875       IF( dzoom.LT.1.)  THEN
    5976         dzoom = dzoom * depi
     
    7390       ENDDO
    7491
    75        DO i = 0, nmax2
     92       DO i = nmax, nmax2
    7693
    7794       fa  = tau*  ( dzoom/2.  - xtild(i) )
    7895       fb  = xtild(i) *  ( pi - xtild(i) )
    7996
    80 
    81        if ( xtild(i) .ne. 0. .and. xtild(i) .ne. pi) then
    8297         IF( 200.* fb .LT. - fa )   THEN
    8398           fhyp ( i) = - 1.
     
    85100           fhyp ( i) =   1.
    86101         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
    88111         ENDIF
    89        endif
    90 
    91        IF ( xtild(i).EQ. 0. )  fhyp(i) =  1.
    92        IF ( xtild(i).EQ. pi )  fhyp(i) = -1.
     112
     113        IF ( xtild(i).EQ. 0. )  fhyp(i) =  1.
     114        IF ( xtild(i).EQ. pi )  fhyp(i) = -1.
    93115
    94116       ENDDO
     
    100122
    101123       DO i = nmax +1,nmax2
     124
     125       xmoy    = 0.5 * ( xtild(i-1) + xtild( i ) )
     126       fa  = tau*  ( dzoom/2.  - xmoy )
     127       fb  = xmoy *  ( pi - xmoy )
     128
     129       IF( 200.* fb .LT. - fa )   THEN
     130         fxm = - 1.
     131       ELSEIF( 200. * fb .LT. fa ) THEN
     132         fxm =   1.
     133       ELSE
     134            IF( ABS(fa).LT.1.e-13.AND.ABS(fb).LT.1.e-13)  THEN
     135                IF(   200.*fb + fa.LT.1.e-10 )  THEN
     136                    fxm   = - 1.
     137                ELSEIF( 200.*fb - fa.LT.1.e-10 )  THEN
     138                    fxm   =   1.
     139                ENDIF
     140            ELSE
     141                    fxm   =  TANH ( fa/fb )
     142            ENDIF
     143       ENDIF
     144
     145       IF ( xmoy.EQ. 0. )  fxm  =  1.
     146       IF ( xmoy.EQ. pi )  fxm  = -1.
     147
     148       ffdx = ffdx + fxm * ( xtild(i) - xtild(i-1) )
     149
     150       ENDDO
     151
     152        beta  = ( grossism * ffdx - pi ) / ( ffdx - pi )
     153
     154       IF( 2.*beta - grossism.LE. 0.)  THEN
     155        WRITE(6,*) ' **  Attention ! La valeur beta calculee dans la rou
     156     ,tine fxhyp est mauvaise ! '
     157        WRITE(6,*)'Modifier les valeurs de  grossismx ,tau ou dzoomx ',
     158     , ' et relancer ! ***  '
     159        CALL ABORT
     160       ENDIF
     161c
     162c   .....  calcul  de  Xprimt   .....
     163c
     164       
     165       DO i = nmax, nmax2
     166        Xprimt(i) = beta  + ( grossism - beta ) * fhyp(i)
     167       ENDDO
     168c   
     169       DO i =  nmax+1, nmax2
     170        Xprimt( nmax2 - i ) = Xprimt( i )
     171       ENDDO
     172c
     173
     174c   .....  Calcul  de  Xf     ........
     175
     176       Xf(0) = - pi
     177
     178       DO i =  nmax +1, nmax2
    102179
    103180       xmoy    = 0.5 * ( xtild(i-1) + xtild( i ) )
     
    113190       ENDIF
    114191
    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        ENDDO
    120 
    121         beta  = ( grossism * ffdx - pi ) / ( ffdx - pi )
    122 
    123        IF( 2.*beta - grossism.LE. 0.)  THEN
    124         WRITE(6,*) ' **  Attention ! La valeur beta calculee dans la rou
    125      ,tine fxhyp est mauvaise ! '
    126         WRITE(6,*)'Modifier les valeurs de  grossismx ,tau ou dzoomx ',
    127      , ' et relancer ! ***  '
    128         CALL ABORT
    129        ENDIF
    130 c
    131 c   .....  calcul  de  Xprimt   .....
    132 c
    133        
    134        DO i = nmax, nmax2
    135         Xprimt(i) = beta  + ( grossism - beta ) * fhyp(i)
    136        ENDDO
    137 c   
    138        DO i =  nmax+1, nmax2
    139         Xprimt( nmax2 - i ) = Xprimt( i )
    140        ENDDO
    141 c
    142 
    143 c   .....  Calcul  de  Xf     ........
    144 
    145        Xf(0) = - pi
    146 
    147        DO i =  nmax +1, nmax2
    148 
    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 )   THEN
    154          fxm = - 1.
    155        ELSEIF( 200. * fb .LT. fa ) THEN
    156          fxm =   1.
    157        ELSE
    158          fxm =  TANH ( fa/fb )
    159        ENDIF
    160 
    161192       IF ( xmoy.EQ. 0. )  fxm =  1.
    162193       IF ( xmoy.EQ. pi )  fxm = -1.
     
    185216
    186217       IF( ik.EQ.1 )        THEN
    187          xuv = - 0.25
     218         xuv =  -0.25
    188219       ELSE IF ( ik.EQ.2 )  THEN
    189220         xuv =   0.
    190221       ELSE IF ( ik.EQ.3 )  THEN
    191          xuv =   0.5
     222         xuv =   0.50
    192223       ELSE IF ( ik.EQ.4 )  THEN
    193224         xuv =   0.25
     
    196227      xo1   = 0.
    197228
    198       DO 1500 i = 1, iim
    199 
    200       xlon2 = - pi + (FLOAT(i) + xuv - 0.75) * depi / FLOAT(iim)
     229      ii1=1
     230      ii2=iim
     231      IF(ik.EQ.1.and.grossism.EQ.1.) THEN
     232        ii1 = 2
     233        ii2 = iim+1
     234      ENDIF
     235
     236      DO 1500 i = ii1, ii2
     237
     238      xlon2 = - pi + (FLOAT(i) + xuv - decalx) * depi / FLOAT(iim)
    201239
    202240      Xfi    = xlon2
     
    2482861500   CONTINUE
    249287
     288
     289       IF(ik.EQ.1.and.grossism.EQ.1.)  THEN
     290         xvrai(1)    = xvrai(iip1)-depi
     291         xxprim(1)   = xxprim(iip1)
     292       ENDIF
    250293       DO i = 1 , iim
    251294        xlon(i)     = xvrai(i)
    252295        xprimm(i)   = xxprim(i)
    253296       ENDDO
    254        
     297 
    255298       DO i = 1, iim -1
    256299        IF( xvrai(i+1). LT. xvrai(i) )  THEN
     
    271314       ENDDO
    272315
    273       IF(champmin .GE. - pi .AND. champmax.LE. pi )  THEN
     316      IF(champmin .GE.-pi-0.10.and.champmax.LE.pi+0.10 )  THEN
    274317                GO TO 1600
    275318      ELSE
     
    301344          IF( ik.EQ.1 )  THEN
    302345           DO i = iim,1,-1
    303             IF( xvrai(i).LE. pi ) GO TO 90
     346             IF( xvrai(i).LE. pi ) GO TO 90
    304347           ENDDO
    305348             WRITE(6,*) ' PBS.  2 ! Xvrai plus grand  que   pi ! '
     
    336379c          WRITE(6,18)
    337380c          WRITE(6,68) xvrai
    338 ccc          WRITE(6,*) ' XPRIM k ',ik
    339 ccc          WRITE(6,566)  xprimm
    340 
    341            DO i = 1,iim + 1
     381c          WRITE(6,*) ' XPRIM k ',ik
     382c          WRITE(6,566)  xprimm
     383
     384           DO i = 1,iim +1
    342385             rlonm025(i) = xlon( i )
    343386            xprimm025(i) = xprimm(i)
    344387           ENDDO
    345 
     388           
    346389         ELSE IF( ik.EQ.2 )  THEN
    347390c          WRITE(6,18)
    348391c          WRITE(6,*)  ' XLON aux pts. V   apres ( en  deg. ) '
    349392c          WRITE(6,68) xvrai
    350 ccc          WRITE(6,*) ' XPRIM k ',ik
    351 ccc          WRITE(6,566)  xprimm
     393c          WRITE(6,*) ' XPRIM k ',ik
     394c          WRITE(6,566)  xprimm
    352395
    353396           DO i = 1,iim + 1
     
    360403c          WRITE(6,*)  ' XLON aux pts. U   apres ( en  deg. ) '
    361404c          WRITE(6,68) xvrai
    362 ccc          WRITE(6,*) ' XPRIM ik ',ik
    363 ccc          WRITE(6,566)  xprimm
     405c          WRITE(6,*) ' XPRIM ik ',ik
     406c          WRITE(6,566)  xprimm
    364407
    365408           DO i = 1,iim + 1
     
    372415c          WRITE(6,*)  ' XLON aux pts. V+0.25   apres ( en  deg. ) '
    373416c          WRITE(6,68) xvrai
    374 ccc          WRITE(6,*) ' XPRIM ik ',ik
    375 ccc          WRITE(6,566)  xprimm
     417c          WRITE(6,*) ' XPRIM ik ',ik
     418c          WRITE(6,566)  xprimm
    376419
    377420           DO i = 1,iim + 1
     
    40344624     FORMAT(2x,'Parametres xzoom,gross,tau ,dzoom pour fxhyp ',4f8.3)
    40444768     FORMAT(1x,7f9.2)
     448566    FORMAT(1x,7f9.4)
    405449
    406450       RETURN
Note: See TracChangeset for help on using the changeset viewer.