source: LMDZ5/trunk/libf/dyn3dpar/fxhyp.F @ 5332

Last change on this file since 5332 was 1939, checked in by lguez, 11 years ago

Same as revision 1930: replaced abort by abort_gcm.

Also replaced real*8 by real(kind=8) (was done way back in revision 1220 for
dyn3d/fxhyp.F but not dyn3dpar/fxhyp.F).

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