source: LMDZ.3.3/branches/rel-1-0-patch/libf/dyn3d/fxhyp.F @ 5230

Last change on this file since 5230 was 298, checked in by lmdz, 23 years ago

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

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 11.7 KB
RevLine 
[2]1       SUBROUTINE fxhyp ( xzoomdeg,grossism,dzoom,tau ,
[212]2     , rlonm025,xprimm025,rlonv,xprimv,rlonu,xprimu,rlonp025,xprimp025,
3     , champmin,champmax                                               )
[2]4
5c      Auteur :  P. Le Van
6
7       IMPLICIT NONE
8
9c    Calcule les longitudes et derivees dans la grille du GCM pour une
10c     fonction f(x) a tangente  hyperbolique  .
11c
12c     grossism etant le grossissement ( = 2 si 2 fois, = 3 si 3 fois,etc.)
13c     dzoom  etant  la distance totale de la zone du zoom
[203]14c     tau  la raideur de la transition de l'interieur a l'exterieur du zoom
[2]15c
[203]16c    On doit avoir grossism x dzoom <  pi ( radians )   , en longitude.
17c   ********************************************************************
[2]18
[203]19
[2]20       INTEGER nmax, nmax2
[203]21       PARAMETER (  nmax = 30000, nmax2 = 2*nmax )
[298]22c
23       LOGICAL scal180
24       PARAMETER ( scal180 = .TRUE. )
[2]25
[298]26c      scal180 = .TRUE.  si on veut avoir le premier point scalaire pour   
27c      une grille reguliere ( grossism = 1.,tau=0.,clon=0. ) a -180. degres.
28c      sinon scal180 = .FALSE.
29
[2]30#include "dimensions.h"
31#include "paramet.h"
32       
33c     ......  arguments  d'entree   .......
34c
35       REAL xzoomdeg,dzoom,tau,grossism
[203]36
37c    ......   arguments  de  sortie  ......
38
[2]39       REAL rlonm025(iip1),xprimm025(iip1),rlonv(iip1),xprimv(iip1),
40     ,  rlonu(iip1),xprimu(iip1),rlonp025(iip1),xprimp025(iip1)
41
[203]42c     .... variables locales  ....
[2]43c
[203]44       REAL*8 xlon(iip1),xprimm(iip1),xuv
45       REAL*8 xtild(0:nmax2)
[216]46       REAL*8 fhyp(0:nmax2),ffdx,beta,Xprimt(0:nmax2)
[203]47       REAL*8 Xf(0:nmax2),xxpr(0:nmax2)
48       REAL*8 xvrai(iip1),xxprim(iip1)
49       REAL*8 pi,depi,epsilon,xzoom,fa,fb
50       REAL*8 Xf1, Xfi , a0,a1,a2,a3,xi2
[298]51       INTEGER i,it,ik,iter,ii,idif,ii1,ii2
[203]52       REAL*8 xi,xo1,xmoy,xlon2,fxm,Xprimin
[298]53       REAL*8 champmin,champmax,decalx
[2]54       INTEGER is2
55       SAVE is2
[212]56
[203]57       REAL*8 heavyside
58       EXTERNAL coefpoly,heavyside
[2]59
60       pi       = 2. * ASIN(1.)
61       depi     = 2. * pi
[203]62       epsilon  = 1.e-3
[2]63       xzoom    = xzoomdeg * pi/180.
[298]64c
65           decalx   = .75
66       IF( grossism.EQ.1..AND.scal180 )  THEN
67           decalx   = 1.
68       ENDIF
[2]69
[298]70       WRITE(6,*) 'FXHYP scal180,decalx', scal180,decalx
71c
[203]72       IF( dzoom.LT.1.)  THEN
73         dzoom = dzoom * depi
74       ELSEIF( dzoom.LT. 25. ) THEN
75         WRITE(6,*) ' Le param. dzoomy pour fxhyp est trop petit ! L aug
76     ,menter et relancer ! '
77         STOP 1
78       ELSE
79         dzoom = dzoom * pi/180.
80       ENDIF
[2]81
[212]82       WRITE(6,*) ' xzoom( rad.),grossism,tau,dzoom (radians)'
83       WRITE(6,24) xzoom,grossism,tau,dzoom
84
[2]85       DO i = 0, nmax2
[203]86        xtild(i) = - pi + FLOAT(i) * depi /nmax2
[2]87       ENDDO
88
[298]89       DO i = nmax, nmax2
[203]90
91       fa  = tau*  ( dzoom/2.  - xtild(i) )
92       fb  = xtild(i) *  ( pi - xtild(i) )
93
[251]94         IF( 200.* fb .LT. - fa )   THEN
95           fhyp ( i) = - 1.
96         ELSEIF( 200. * fb .LT. fa ) THEN
97           fhyp ( i) =   1.
98         ELSE
[298]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
[251]108         ENDIF
[203]109
[298]110        IF ( xtild(i).EQ. 0. )  fhyp(i) =  1.
111        IF ( xtild(i).EQ. pi )  fhyp(i) = -1.
[203]112
[2]113       ENDDO
114
115cc  ....  Calcul  de  beta  ....
[203]116c   ............................
[2]117
[203]118       ffdx = 0.
119
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
[298]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
[203]140       ENDIF
141
[298]142       IF ( xmoy.EQ. 0. )  fxm  =  1.
143       IF ( xmoy.EQ. pi )  fxm  = -1.
144
[203]145       ffdx = ffdx + fxm * ( xtild(i) - xtild(i-1) )
146
[2]147       ENDDO
148
[203]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
[2]158c
159c   .....  calcul  de  Xprimt   .....
160c
161       
[203]162       DO i = nmax, nmax2
[2]163        Xprimt(i) = beta  + ( grossism - beta ) * fhyp(i)
164       ENDDO
165c   
[203]166       DO i =  nmax+1, nmax2
[2]167        Xprimt( nmax2 - i ) = Xprimt( i )
168       ENDDO
169c
170
171c   .....  Calcul  de  Xf     ........
172
[203]173       Xf(0) = - pi
174
175       DO i =  nmax +1, nmax2
176
177       xmoy    = 0.5 * ( xtild(i-1) + xtild( i ) )
178       fa  = tau*  ( dzoom/2.  - xmoy )
179       fb  = xmoy *  ( pi - xmoy )
180
181       IF( 200.* fb .LT. - fa )   THEN
182         fxm = - 1.
183       ELSEIF( 200. * fb .LT. fa ) THEN
184         fxm =   1.
185       ELSE
186         fxm =  TANH ( fa/fb )
187       ENDIF
188
189       IF ( xmoy.EQ. 0. )  fxm =  1.
190       IF ( xmoy.EQ. pi )  fxm = -1.
191       xxpr(i)    = beta + ( grossism - beta ) * fxm
192
[2]193       ENDDO
194
[203]195       DO i = nmax+1, nmax2
196        xxpr(nmax2-i+1) = xxpr(i)
[2]197       ENDDO
198
199        DO i=1,nmax2
200         Xf(i)   = Xf(i-1) + xxpr(i) * ( xtild(i) - xtild(i-1) )
201        ENDDO
202
203
204c    *****************************************************************
205c
206
207c     .....  xuv = 0.   si  calcul  aux pts   scalaires   ........
208c     .....  xuv = 0.5  si  calcul  aux pts      U        ........
209c
[212]210      WRITE(6,18)
[2]211c
212      DO 5000  ik = 1, 4
213
214       IF( ik.EQ.1 )        THEN
[298]215         xuv =  -0.25
[2]216       ELSE IF ( ik.EQ.2 )  THEN
217         xuv =   0.
218       ELSE IF ( ik.EQ.3 )  THEN
[298]219         xuv =   0.50
[2]220       ELSE IF ( ik.EQ.4 )  THEN
221         xuv =   0.25
222       ENDIF
223
[203]224      xo1   = 0.
[2]225
[298]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
[2]232
[298]233      DO 1500 i = ii1, ii2
[2]234
[298]235      xlon2 = - pi + (FLOAT(i) + xuv - decalx) * depi / FLOAT(iim)
236
[203]237      Xfi    = xlon2
[2]238c
239      DO 250 it =  nmax2,0,-1
[203]240      IF( Xfi.GE.Xf(it))  GO TO 350
[2]241250   CONTINUE
242
243      it = 0
244
245350   CONTINUE
246
[203]247c    ......  Calcul de   Xf(xi)    ......
248c
249      xi  = xtild(it)
[2]250
[203]251      IF(it.EQ.nmax2)  THEN
252       it       = nmax2 -1
253       Xf(it+1) = pi
254      ENDIF
255c  .....................................................................
256c
257c   Appel de la routine qui calcule les coefficients a0,a1,a2,a3 d'un
258c   polynome de degre 3  qui passe  par les points (Xf(it),xtild(it) )
259c          et (Xf(it+1),xtild(it+1) )
[2]260
[203]261       CALL coefpoly ( Xf(it),Xf(it+1),Xprimt(it),Xprimt(it+1),
262     ,                xtild(it),xtild(it+1),  a0, a1, a2, a3  )
263
264       Xf1     = Xf(it)
265       Xprimin = a1 + 2.* a2 * xi + 3.*a3 * xi *xi
266
267       DO 500 iter = 1,300
268        xi = xi - ( Xf1 - Xfi )/ Xprimin
269
270        IF( ABS(xi-xo1).LE.epsilon)  GO TO 550
271         xo1      = xi
272         xi2      = xi * xi
273         Xf1      = a0 +  a1 * xi +     a2 * xi2  +     a3 * xi2 * xi
274         Xprimin  =       a1      + 2.* a2 *  xi  + 3.* a3 * xi2
[2]275500   CONTINUE
[203]276        WRITE(6,*) ' Pas de solution ***** ',i,xlon2,iter
277          STOP 6
[2]278550   CONTINUE
279
[203]280       xxprim(i) = depi/ ( FLOAT(iim) * Xprimin )
281       xvrai(i)  =  xi + xzoom
[2]282
2831500   CONTINUE
284
[298]285
286       IF(ik.EQ.1.and.grossism.EQ.1.)  THEN
287         xvrai(1)    = xvrai(iip1)-depi
288         xxprim(1)   = xxprim(iip1)
289       ENDIF
[2]290       DO i = 1 , iim
[203]291        xlon(i)     = xvrai(i)
[2]292        xprimm(i)   = xxprim(i)
293       ENDDO
[298]294 
[2]295       DO i = 1, iim -1
296        IF( xvrai(i+1). LT. xvrai(i) )  THEN
[203]297         WRITE(6,*) ' PBS. avec rlonu(',i+1,') plus petit que rlonu(',i,
298     ,  ')'
299        STOP 7
[2]300        ENDIF
301       ENDDO
302c
303c   ... Reorganisation  des  longitudes  pour les avoir  entre - pi et pi ..
304c   ........................................................................
305
306       champmin =  1.e12
307       champmax = -1.e12
308       DO i = 1, iim
309        champmin = MIN( champmin,xvrai(i) )
310        champmax = MAX( champmax,xvrai(i) )
311       ENDDO
312
[298]313      IF(champmin .GE.-pi-0.10.and.champmax.LE.pi+0.10 )  THEN
[2]314                GO TO 1600
315      ELSE
[203]316       WRITE(6,*) 'Reorganisation des longitudes pour avoir entre - pi',
317     ,  ' et pi '
[2]318c
319        IF( xzoom.LE.0.)  THEN
320          IF( ik.EQ. 1 )  THEN
321          DO i = 1, iim
322           IF( xvrai(i).GE. - pi )  GO TO 80
323          ENDDO
[203]324            WRITE(6,*)  ' PBS. 1 !  Xvrai plus petit que  - pi ! '
325            STOP 8
[2]326 80       CONTINUE
327          is2 = i
328          ENDIF
329
330          IF( is2.NE. 1 )  THEN
331            DO ii = is2 , iim
332             xlon  (ii-is2+1) = xvrai(ii)
333             xprimm(ii-is2+1) = xxprim(ii)
334            ENDDO
335            DO ii = 1 , is2 -1
336             xlon  (ii+iim-is2+1) = xvrai(ii) + depi
337             xprimm(ii+iim-is2+1) = xxprim(ii)
338            ENDDO
339          ENDIF
340        ELSE
341          IF( ik.EQ.1 )  THEN
342           DO i = iim,1,-1
[298]343             IF( xvrai(i).LE. pi ) GO TO 90
[2]344           ENDDO
[203]345             WRITE(6,*) ' PBS.  2 ! Xvrai plus grand  que   pi ! '
346              STOP 9
[2]347 90        CONTINUE
348            is2 = i
349          ENDIF
350           idif = iim -is2
351           DO ii = 1, is2
352            xlon  (ii+idif) = xvrai(ii)
353            xprimm(ii+idif) = xxprim(ii)
354           ENDDO
355           DO ii = 1, idif
356            xlon (ii)  = xvrai (ii+is2) - depi
357            xprimm(ii) = xxprim(ii+is2)
358           ENDDO
359         ENDIF
360      ENDIF
361c
362c     .........   Fin  de la reorganisation   ............................
363
364 1600    CONTINUE
365
366
367         xlon  ( iip1)  = xlon(1) + depi
368         xprimm( iip1 ) = xprimm (1 )
369       
370         DO i = 1, iim+1
371         xvrai(i) = xlon(i)*180./pi
372         ENDDO
373
[203]374         IF( ik.EQ.1 )  THEN
[212]375c          WRITE(6,*)  ' XLON aux pts. V-0.25   apres ( en  deg. ) '
376c          WRITE(6,18)
377c          WRITE(6,68) xvrai
[298]378c          WRITE(6,*) ' XPRIM k ',ik
379c          WRITE(6,566)  xprimm
[2]380
[298]381           DO i = 1,iim +1
[2]382             rlonm025(i) = xlon( i )
383            xprimm025(i) = xprimm(i)
384           ENDDO
[298]385           
[2]386         ELSE IF( ik.EQ.2 )  THEN
[212]387c          WRITE(6,18)
388c          WRITE(6,*)  ' XLON aux pts. V   apres ( en  deg. ) '
389c          WRITE(6,68) xvrai
[298]390c          WRITE(6,*) ' XPRIM k ',ik
391c          WRITE(6,566)  xprimm
[203]392
[2]393           DO i = 1,iim + 1
394             rlonv(i) = xlon( i )
395            xprimv(i) = xprimm(i)
396           ENDDO
[203]397
398         ELSE IF( ik.EQ.3)   THEN
[212]399c          WRITE(6,18)
400c          WRITE(6,*)  ' XLON aux pts. U   apres ( en  deg. ) '
401c          WRITE(6,68) xvrai
[298]402c          WRITE(6,*) ' XPRIM ik ',ik
403c          WRITE(6,566)  xprimm
[203]404
[2]405           DO i = 1,iim + 1
406             rlonu(i) = xlon( i )
407            xprimu(i) = xprimm(i)
408           ENDDO
[203]409
[2]410         ELSE IF( ik.EQ.4 )  THEN
[212]411c          WRITE(6,18)
412c          WRITE(6,*)  ' XLON aux pts. V+0.25   apres ( en  deg. ) '
413c          WRITE(6,68) xvrai
[298]414c          WRITE(6,*) ' XPRIM ik ',ik
415c          WRITE(6,566)  xprimm
[203]416
[2]417           DO i = 1,iim + 1
418             rlonp025(i) = xlon( i )
419            xprimp025(i) = xprimm(i)
420           ENDDO
[203]421
[2]422         ENDIF
423
4245000    CONTINUE
425c
[212]426       WRITE(6,18)
427c
[203]428c    ...........  fin  de la boucle  do 5000      ............
[2]429
[203]430        DO i = 1, iim
431         xlon(i) = rlonv(i+1) - rlonv(i)
[2]432        ENDDO
[203]433        champmin =  1.e12
434        champmax = -1.e12
[2]435        DO i = 1, iim
[203]436         champmin = MIN( champmin, xlon(i) )
437         champmax = MAX( champmax, xlon(i) )
[2]438        ENDDO
[203]439         champmin = champmin * 180./pi
440         champmax = champmax * 180./pi
[2]441
[203]44218     FORMAT(/)
[212]44324     FORMAT(2x,'Parametres xzoom,gross,tau ,dzoom pour fxhyp ',4f8.3)
[203]44468     FORMAT(1x,7f9.2)
[298]445566    FORMAT(1x,7f9.4)
[2]446
[203]447       RETURN
448       END
Note: See TracBrowser for help on using the repository browser.