source: LMDZ.3.3/trunk/libf/dyn3d/fxhyp.F @ 346

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