/[lmdze]/trunk/dyn3d/fxhyp.f
ViewVC logotype

Diff of /trunk/dyn3d/fxhyp.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/dyn3d/fxhyp.f revision 76 by guez, Fri Nov 15 18:45:49 2013 UTC trunk/Sources/dyn3d/fxhyp.f revision 151 by guez, Tue Jun 23 15:14:20 2015 UTC
# Line 1  Line 1 
1  !  module fxhyp_m
 ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/fxhyp.F,v 1.2 2005/06/03 09:11:32 fairhead Exp $  
 !  
 c  
 c  
        SUBROUTINE fxhyp ( xzoomdeg,grossism,dzooma,tau ,  
      , rlonm025,xprimm025,rlonv,xprimv,rlonu,xprimu,rlonp025,xprimp025,  
      , champmin,champmax                                               )  
   
 c      Auteur :  P. Le Van  
   
        use dimens_m  
       use paramet_m  
        IMPLICIT NONE  
   
 c    Calcule les longitudes et derivees dans la grille du GCM pour une  
 c     fonction f(x) a tangente  hyperbolique  .  
 c  
 c     grossism etant le grossissement ( = 2 si 2 fois, = 3 si 3 fois,etc.)  
 c     dzoom  etant  la distance totale de la zone du zoom  
 c     tau  la raideur de la transition de l'interieur a l'exterieur du zoom  
 c  
 c    On doit avoir grossism x dzoom <  pi ( radians )   , en longitude.  
 c   ********************************************************************  
   
   
        INTEGER nmax, nmax2  
        PARAMETER (  nmax = 30000, nmax2 = 2*nmax )  
 c  
        LOGICAL scal180  
        PARAMETER ( scal180 = .TRUE. )  
   
 c      scal180 = .TRUE.  si on veut avoir le premier point scalaire pour    
 c      une grille reguliere ( grossism = 1.,tau=0.,clon=0. ) a -180. degres.  
 c      sinon scal180 = .FALSE.  
   
         
 c     ......  arguments  d'entree   .......  
 c  
        REAL xzoomdeg,dzooma,tau,grossism  
   
 c    ......   arguments  de  sortie  ......  
   
        REAL rlonm025(iip1),xprimm025(iip1),rlonv(iip1),xprimv(iip1),  
      ,  rlonu(iip1),xprimu(iip1),rlonp025(iip1),xprimp025(iip1)  
   
 c     .... variables locales  ....  
 c  
        REAL   dzoom  
        DOUBLE PRECISION xlon(iip1),xprimm(iip1),xuv  
        DOUBLE PRECISION xtild(0:nmax2)  
        DOUBLE PRECISION fhyp(0:nmax2),ffdx,beta,Xprimt(0:nmax2)  
        DOUBLE PRECISION Xf(0:nmax2),xxpr(0:nmax2)  
        DOUBLE PRECISION xvrai(iip1),xxprim(iip1)  
        DOUBLE PRECISION pi,depi,epsilon,xzoom,fa,fb  
        DOUBLE PRECISION Xf1, Xfi , a0,a1,a2,a3,xi2  
        INTEGER i,it,ik,iter,ii,idif,ii1,ii2  
        DOUBLE PRECISION xi,xo1,xmoy,xlon2,fxm,Xprimin  
        DOUBLE PRECISION champmin,champmax,decalx  
        INTEGER is2  
        SAVE is2  
   
        DOUBLE PRECISION heavyside  
   
        pi       = 2. * ASIN(1.)  
        depi     = 2. * pi  
        epsilon  = 1.e-3  
        xzoom    = xzoomdeg * pi/180.  
 c  
            decalx   = .75  
        IF( grossism.EQ.1..AND.scal180 )  THEN  
            decalx   = 1.  
        ENDIF  
   
        WRITE(6,*) 'FXHYP scal180,decalx', scal180,decalx  
 c  
        IF( dzooma.LT.1.)  THEN  
          dzoom = dzooma * depi  
        ELSEIF( dzooma.LT. 25. ) THEN  
          WRITE(6,*) ' Le param. dzoomx pour fxhyp est trop petit ! L aug  
      ,menter et relancer ! '  
          STOP 1  
        ELSE  
          dzoom = dzooma * pi/180.  
        ENDIF  
2    
3         WRITE(6,*) ' xzoom( rad.),grossism,tau,dzoom (radians)'    IMPLICIT NONE
        WRITE(6,24) xzoom,grossism,tau,dzoom  
4    
5         DO i = 0, nmax2  contains
6          xtild(i) = - pi + FLOAT(i) * depi /nmax2  
7         ENDDO    SUBROUTINE fxhyp(xprimm025, rlonv, xprimv, rlonu, xprimu, xprimp025)
8    
9         DO i = nmax, nmax2      ! From LMDZ4/libf/dyn3d/fxhyp.F, version 1.2, 2005/06/03 09:11:32
10        ! Author: P. Le Van, from formulas by R. Sadourny
11         fa  = tau*  ( dzoom/2.  - xtild(i) )  
12         fb  = xtild(i) *  ( pi - xtild(i) )      ! Calcule les longitudes et dérivées dans la grille du GCM pour
13        ! une fonction x_f(\tilde x) à dérivée tangente hyperbolique.
14           IF( 200.* fb .LT. - fa )   THEN  
15             fhyp ( i) = - 1.      ! Il vaut mieux avoir : grossismx \times delta < pi
16           ELSEIF( 200. * fb .LT. fa ) THEN  
17             fhyp ( i) =   1.      ! Le premier point scalaire pour une grille regulière (grossismx =
18           ELSE      ! 1) avec clon = 0 est à - 180 degrés.
19              IF( ABS(fa).LT.1.e-13.AND.ABS(fb).LT.1.e-13)  THEN  
20                  IF(   200.*fb + fa.LT.1.e-10 )  THEN      USE dimens_m, ONLY: iim
21                      fhyp ( i ) = - 1.      use dynetat0_m, only: clon, grossismx, dzoomx, taux
22                  ELSEIF( 200.*fb - fa.LT.1.e-10 )  THEN      use invert_zoom_x_m, only: invert_zoom_x, nmax
23                      fhyp ( i )  =   1.      use nr_util, only: pi, pi_d, twopi, twopi_d, arth
24                  ENDIF      use principal_cshift_m, only: principal_cshift
25              ELSE      use tanh_cautious_m, only: tanh_cautious
26                      fhyp ( i )  =  TANH ( fa/fb )  
27              ENDIF      REAL, intent(out):: xprimm025(:), rlonv(:), xprimv(:) ! (iim + 1)
28           ENDIF      real, intent(out):: rlonu(:), xprimu(:), xprimp025(:) ! (iim + 1)
29          IF ( xtild(i).EQ. 0. )  fhyp(i) =  1.  
30          IF ( xtild(i).EQ. pi )  fhyp(i) = -1.      ! Local:
31        real rlonm025(iim + 1), rlonp025(iim + 1), d_rlonv(iim)
32         ENDDO      REAL delta, step
33        DOUBLE PRECISION, dimension(0:nmax):: xtild, fhyp, G, Xf, ffdx
34  cc  ....  Calcul  de  beta  ....      DOUBLE PRECISION beta
35        INTEGER i, is2
36         ffdx = 0.      DOUBLE PRECISION xmoy(nmax), fxm(nmax)
37    
38         DO i = nmax +1,nmax2      !----------------------------------------------------------------------
39    
40         xmoy    = 0.5 * ( xtild(i-1) + xtild( i ) )      print *, "Call sequence information: fxhyp"
41         fa  = tau*  ( dzoom/2.  - xmoy )  
42         fb  = xmoy *  ( pi - xmoy )      if (grossismx == 1.) then
43           step = twopi / iim
44         IF( 200.* fb .LT. - fa )   THEN  
45           fxm = - 1.         xprimm025(:iim) = step
46         ELSEIF( 200. * fb .LT. fa ) THEN         xprimp025(:iim) = step
47           fxm =   1.         xprimv(:iim) = step
48         ELSE         xprimu(:iim) = step
49              IF( ABS(fa).LT.1.e-13.AND.ABS(fb).LT.1.e-13)  THEN  
50                  IF(   200.*fb + fa.LT.1.e-10 )  THEN         rlonv(:iim) = arth(- pi + clon, step, iim)
51                      fxm   = - 1.         rlonm025(:iim) = rlonv(:iim) - 0.25 * step
52                  ELSEIF( 200.*fb - fa.LT.1.e-10 )  THEN         rlonp025(:iim) = rlonv(:iim) + 0.25 * step
53                      fxm   =   1.         rlonu(:iim) = rlonv(:iim) + 0.5 * step
54                  ENDIF      else
55              ELSE         delta = dzoomx * twopi_d
56                      fxm   =  TANH ( fa/fb )         xtild = arth(0d0, pi_d / nmax, nmax + 1)
57              ENDIF         forall (i = 1:nmax) xmoy(i) = 0.5d0 * (xtild(i-1) + xtild(i))
58         ENDIF  
59           ! Compute fhyp:
60         IF ( xmoy.EQ. 0. )  fxm  =  1.         fhyp(1:nmax - 1) = tanh_cautious(taux * (delta / 2d0 &
61         IF ( xmoy.EQ. pi )  fxm  = -1.              - xtild(1:nmax - 1)), xtild(1:nmax - 1) &
62                * (pi_d - xtild(1:nmax - 1)))
63         ffdx = ffdx + fxm * ( xtild(i) - xtild(i-1) )         fhyp(0) = 1d0
64           fhyp(nmax) = -1d0
65         ENDDO  
66           fxm = tanh_cautious(taux * (delta / 2d0 - xmoy), xmoy * (pi_d - xmoy))
67          beta  = ( grossism * ffdx - pi ) / ( ffdx - pi )  
68           ! Compute \int_0 ^{\tilde x} F:
69         IF( 2.*beta - grossism.LE. 0.)  THEN  
70          WRITE(6,*) ' **  Attention ! La valeur beta calculee dans la rou         ffdx(0) = 0d0
71       ,tine fxhyp est mauvaise ! '  
72          WRITE(6,*)'Modifier les valeurs de  grossismx ,tau ou dzoomx ',         DO i = 1, nmax
73       , ' et relancer ! ***  '            ffdx(i) = ffdx(i - 1) + fxm(i) * (xtild(i) - xtild(i-1))
74          STOP 1         END DO
75         ENDIF  
76  c         print *, "ffdx(nmax) = ", ffdx(nmax)
77  c   .....  calcul  de  Xprimt   .....         beta = (pi_d - grossismx * ffdx(nmax)) / (pi_d - ffdx(nmax))
78  c         print *, "beta = ", beta
79          
80         DO i = nmax, nmax2         IF (2d0 * beta - grossismx <= 0d0) THEN
81          Xprimt(i) = beta  + ( grossism - beta ) * fhyp(i)            print *, 'Bad choice of grossismx, taux, dzoomx.'
82         ENDDO            print *, 'Decrease dzoomx or grossismx.'
83  c              STOP 1
84         DO i =  nmax+1, nmax2         END IF
85          Xprimt( nmax2 - i ) = Xprimt( i )  
86         ENDDO         G = beta + (grossismx - beta) * fhyp
87  c  
88           Xf(:nmax - 1) = beta * xtild(:nmax - 1) + (grossismx - beta) &
89  c   .....  Calcul  de  Xf     ........              * ffdx(:nmax - 1)
90           Xf(nmax) = pi_d
91         Xf(0) = - pi  
92           call invert_zoom_x(xf, xtild, G, rlonm025(:iim), xprimm025(:iim), &
93         DO i =  nmax +1, nmax2              xuv = - 0.25d0)
94           call invert_zoom_x(xf, xtild, G, rlonv(:iim), xprimv(:iim), xuv = 0d0)
95         xmoy    = 0.5 * ( xtild(i-1) + xtild( i ) )         call invert_zoom_x(xf, xtild, G, rlonu(:iim), xprimu(:iim), xuv = 0.5d0)
96         fa  = tau*  ( dzoom/2.  - xmoy )         call invert_zoom_x(xf, xtild, G, rlonp025(:iim), xprimp025(:iim), &
97         fb  = xmoy *  ( pi - xmoy )              xuv = 0.25d0)
98        end if
99         IF( 200.* fb .LT. - fa )   THEN  
100           fxm = - 1.      is2 = 0
101         ELSEIF( 200. * fb .LT. fa ) THEN  
102           fxm =   1.      IF (MINval(rlonm025(:iim)) < - pi - 0.1 &
103             .or. MAXval(rlonm025(:iim)) > pi + 0.1) THEN
104           IF (clon <= 0.) THEN
105              is2 = 1
106    
107              do while (rlonm025(is2) < - pi .and. is2 < iim)
108                 is2 = is2 + 1
109              end do
110    
111              if (rlonm025(is2) < - pi) then
112                 print *, 'Rlonm025 plus petit que - pi !'
113                 STOP 1
114              end if
115         ELSE         ELSE
116           fxm =  TANH ( fa/fb )            is2 = iim
117         ENDIF  
118              do while (rlonm025(is2) > pi .and. is2 > 1)
119                 is2 = is2 - 1
120              end do
121    
122              if (rlonm025(is2) > pi) then
123                 print *, 'Rlonm025 plus grand que pi !'
124                 STOP 1
125              end if
126           END IF
127        END IF
128    
129        call principal_cshift(is2, rlonm025, xprimm025)
130        call principal_cshift(is2, rlonv, xprimv)
131        call principal_cshift(is2, rlonu, xprimu)
132        call principal_cshift(is2, rlonp025, xprimp025)
133    
134        forall (i = 1: iim) d_rlonv(i) = rlonv(i + 1) - rlonv(i)
135        print *, "Minimum longitude step:", MINval(d_rlonv) * 180. / pi, "degrees"
136        print *, "Maximum longitude step:", MAXval(d_rlonv) * 180. / pi, "degrees"
137    
138        ! Check that rlonm025 <= rlonv <= rlonp025 <= rlonu:
139        DO i = 1, iim + 1
140           IF (rlonp025(i) < rlonv(i)) THEN
141              print *, 'rlonp025(', i, ') = ', rlonp025(i)
142              print *, "< rlonv(", i, ") = ", rlonv(i)
143              STOP 1
144           END IF
145    
146           IF (rlonv(i) < rlonm025(i)) THEN
147              print *, 'rlonv(', i, ') = ', rlonv(i)
148              print *, "< rlonm025(", i, ") = ", rlonm025(i)
149              STOP 1
150           END IF
151    
152           IF (rlonp025(i) > rlonu(i)) THEN
153              print *, 'rlonp025(', i, ') = ', rlonp025(i)
154              print *, "> rlonu(", i, ") = ", rlonu(i)
155              STOP 1
156           END IF
157        END DO
158    
159         IF ( xmoy.EQ. 0. )  fxm =  1.    END SUBROUTINE fxhyp
        IF ( xmoy.EQ. pi )  fxm = -1.  
        xxpr(i)    = beta + ( grossism - beta ) * fxm  
   
        ENDDO  
   
        DO i = nmax+1, nmax2  
         xxpr(nmax2-i+1) = xxpr(i)  
        ENDDO  
   
         DO i=1,nmax2  
          Xf(i)   = Xf(i-1) + xxpr(i) * ( xtild(i) - xtild(i-1) )  
         ENDDO  
   
   
 c    *****************************************************************  
 c  
   
 c     .....  xuv = 0.   si  calcul  aux pts   scalaires   ........  
 c     .....  xuv = 0.5  si  calcul  aux pts      U        ........  
 c  
       WRITE(6,18)  
 c  
       DO 5000  ik = 1, 4  
   
        IF( ik.EQ.1 )        THEN  
          xuv =  -0.25  
        ELSE IF ( ik.EQ.2 )  THEN  
          xuv =   0.  
        ELSE IF ( ik.EQ.3 )  THEN  
          xuv =   0.50  
        ELSE IF ( ik.EQ.4 )  THEN  
          xuv =   0.25  
        ENDIF  
   
       xo1   = 0.  
   
       ii1=1  
       ii2=iim  
       IF(ik.EQ.1.and.grossism.EQ.1.) THEN  
         ii1 = 2  
         ii2 = iim+1  
       ENDIF  
       DO 1500 i = ii1, ii2  
   
       xlon2 = - pi + (FLOAT(i) + xuv - decalx) * depi / FLOAT(iim)  
   
       Xfi    = xlon2  
 c  
       DO 250 it =  nmax2,0,-1  
       IF( Xfi.GE.Xf(it))  GO TO 350  
 250   CONTINUE  
   
       it = 0  
   
 350   CONTINUE  
   
 c    ......  Calcul de   Xf(xi)    ......  
 c  
       xi  = xtild(it)  
   
       IF(it.EQ.nmax2)  THEN  
        it       = nmax2 -1  
        Xf(it+1) = pi  
       ENDIF  
 c  .....................................................................  
 c  
 c   Appel de la routine qui calcule les coefficients a0,a1,a2,a3 d'un  
 c   polynome de degre 3  qui passe  par les points (Xf(it),xtild(it) )  
 c          et (Xf(it+1),xtild(it+1) )  
   
        CALL coefpoly ( Xf(it),Xf(it+1),Xprimt(it),Xprimt(it+1),  
      ,                xtild(it),xtild(it+1),  a0, a1, a2, a3  )  
   
        Xf1     = Xf(it)  
        Xprimin = a1 + 2.* a2 * xi + 3.*a3 * xi *xi  
   
        DO 500 iter = 1,300  
         xi = xi - ( Xf1 - Xfi )/ Xprimin  
   
         IF( ABS(xi-xo1).LE.epsilon)  GO TO 550  
          xo1      = xi  
          xi2      = xi * xi  
          Xf1      = a0 +  a1 * xi +     a2 * xi2  +     a3 * xi2 * xi  
          Xprimin  =       a1      + 2.* a2 *  xi  + 3.* a3 * xi2  
 500   CONTINUE  
         WRITE(6,*) ' Pas de solution ***** ',i,xlon2,iter  
           STOP 6  
 550   CONTINUE  
   
        xxprim(i) = depi/ ( FLOAT(iim) * Xprimin )  
        xvrai(i)  =  xi + xzoom  
   
 1500   CONTINUE  
   
   
        IF(ik.EQ.1.and.grossism.EQ.1.)  THEN  
          xvrai(1)    = xvrai(iip1)-depi  
          xxprim(1)   = xxprim(iip1)  
        ENDIF  
        DO i = 1 , iim  
         xlon(i)     = xvrai(i)  
         xprimm(i)   = xxprim(i)  
        ENDDO  
        DO i = 1, iim -1  
         IF( xvrai(i+1). LT. xvrai(i) )  THEN  
          WRITE(6,*) ' PBS. avec rlonu(',i+1,') plus petit que rlonu(',i,  
      ,  ')'  
         STOP 7  
         ENDIF  
        ENDDO  
 c  
 c   ... Reorganisation  des  longitudes  pour les avoir  entre - pi et pi ..  
 c   ........................................................................  
   
        champmin =  1.e12  
        champmax = -1.e12  
        DO i = 1, iim  
         champmin = MIN( champmin,xvrai(i) )  
         champmax = MAX( champmax,xvrai(i) )  
        ENDDO  
   
       IF(champmin .GE.-pi-0.10.and.champmax.LE.pi+0.10 )  THEN  
                 GO TO 1600  
       ELSE  
        WRITE(6,*) 'Reorganisation des longitudes pour avoir entre - pi',  
      ,  ' et pi '  
 c  
         IF( xzoom.LE.0.)  THEN  
           IF( ik.EQ. 1 )  THEN  
           DO i = 1, iim  
            IF( xvrai(i).GE. - pi )  GO TO 80  
           ENDDO  
             WRITE(6,*)  ' PBS. 1 !  Xvrai plus petit que  - pi ! '  
             STOP 8  
  80       CONTINUE  
           is2 = i  
           ENDIF  
   
           IF( is2.NE. 1 )  THEN  
             DO ii = is2 , iim  
              xlon  (ii-is2+1) = xvrai(ii)  
              xprimm(ii-is2+1) = xxprim(ii)  
             ENDDO  
             DO ii = 1 , is2 -1  
              xlon  (ii+iim-is2+1) = xvrai(ii) + depi  
              xprimm(ii+iim-is2+1) = xxprim(ii)  
             ENDDO  
           ENDIF  
         ELSE  
           IF( ik.EQ.1 )  THEN  
            DO i = iim,1,-1  
              IF( xvrai(i).LE. pi ) GO TO 90  
            ENDDO  
              WRITE(6,*) ' PBS.  2 ! Xvrai plus grand  que   pi ! '  
               STOP 9  
  90        CONTINUE  
             is2 = i  
           ENDIF  
            idif = iim -is2  
            DO ii = 1, is2  
             xlon  (ii+idif) = xvrai(ii)  
             xprimm(ii+idif) = xxprim(ii)  
            ENDDO  
            DO ii = 1, idif  
             xlon (ii)  = xvrai (ii+is2) - depi  
             xprimm(ii) = xxprim(ii+is2)  
            ENDDO  
          ENDIF  
       ENDIF  
 c  
 c     .........   Fin  de la reorganisation   ............................  
   
  1600    CONTINUE  
   
   
          xlon  ( iip1)  = xlon(1) + depi  
          xprimm( iip1 ) = xprimm (1 )  
         
          DO i = 1, iim+1  
          xvrai(i) = xlon(i)*180./pi  
          ENDDO  
   
          IF( ik.EQ.1 )  THEN  
 c          WRITE(6,*)  ' XLON aux pts. V-0.25   apres ( en  deg. ) '  
 c          WRITE(6,18)  
 c          WRITE(6,68) xvrai  
 c          WRITE(6,*) ' XPRIM k ',ik  
 c          WRITE(6,566)  xprimm  
   
            DO i = 1,iim +1  
              rlonm025(i) = xlon( i )  
             xprimm025(i) = xprimm(i)  
            ENDDO  
          ELSE IF( ik.EQ.2 )  THEN  
 c          WRITE(6,18)  
 c          WRITE(6,*)  ' XLON aux pts. V   apres ( en  deg. ) '  
 c          WRITE(6,68) xvrai  
 c          WRITE(6,*) ' XPRIM k ',ik  
 c          WRITE(6,566)  xprimm  
   
            DO i = 1,iim + 1  
              rlonv(i) = xlon( i )  
             xprimv(i) = xprimm(i)  
            ENDDO  
   
          ELSE IF( ik.EQ.3)   THEN  
 c          WRITE(6,18)  
 c          WRITE(6,*)  ' XLON aux pts. U   apres ( en  deg. ) '  
 c          WRITE(6,68) xvrai  
 c          WRITE(6,*) ' XPRIM ik ',ik  
 c          WRITE(6,566)  xprimm  
   
            DO i = 1,iim + 1  
              rlonu(i) = xlon( i )  
             xprimu(i) = xprimm(i)  
            ENDDO  
   
          ELSE IF( ik.EQ.4 )  THEN  
 c          WRITE(6,18)  
 c          WRITE(6,*)  ' XLON aux pts. V+0.25   apres ( en  deg. ) '  
 c          WRITE(6,68) xvrai  
 c          WRITE(6,*) ' XPRIM ik ',ik  
 c          WRITE(6,566)  xprimm  
   
            DO i = 1,iim + 1  
              rlonp025(i) = xlon( i )  
             xprimp025(i) = xprimm(i)  
            ENDDO  
   
          ENDIF  
   
 5000    CONTINUE  
 c  
        WRITE(6,18)  
 c  
 c    ...........  fin  de la boucle  do 5000      ............  
   
         DO i = 1, iim  
          xlon(i) = rlonv(i+1) - rlonv(i)  
         ENDDO  
         champmin =  1.e12  
         champmax = -1.e12  
         DO i = 1, iim  
          champmin = MIN( champmin, xlon(i) )  
          champmax = MAX( champmax, xlon(i) )  
         ENDDO  
          champmin = champmin * 180./pi  
          champmax = champmax * 180./pi  
   
 18     FORMAT(/)  
 24     FORMAT(2x,'Parametres xzoom,gross,tau ,dzoom pour fxhyp ',4f8.3)  
 68     FORMAT(1x,7f9.2)  
 566    FORMAT(1x,7f9.4)  
160    
161         RETURN  end module fxhyp_m
        END  

Legend:
Removed from v.76  
changed lines
  Added in v.151

  ViewVC Help
Powered by ViewVC 1.1.21