source: LMDZ5/trunk/libf/dyn3d/fxhyp.F @ 1934

Last change on this file since 1934 was 1930, checked in by lguez, 11 years ago

abort, dfloat and pause are not in the Fortran standard. Replaced
abort by abort_gcm and dfloat by dble. Note: I modified dyn3dpar files
that were identical to dyn3d modified files.

  • 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
Line 
1!
2! $Id: fxhyp.F 1930 2014-01-17 16:45:09Z jghattas $
3!
4c
5c
6       SUBROUTINE fxhyp ( xzoomdeg,grossism,dzooma,tau ,
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
40       REAL xzoomdeg,dzooma,tau,grossism
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
49       REAL   dzoom
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
57       INTEGER i,it,ik,iter,ii,idif,ii1,ii2
58       REAL(KIND=8) xi,xo1,xmoy,xlon2,fxm,Xprimin
59       REAL(KIND=8) champmin,champmax,decalx
60       INTEGER is2
61       SAVE is2
62
63       REAL(KIND=8) heavyside
64
65       pi       = 2. * ASIN(1.)
66       depi     = 2. * pi
67       epsilon  = 1.e-3
68       xzoom    = xzoomdeg * pi/180.
69c
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
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
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
102     ,menter et relancer ! '
103         STOP 1
104       ELSE
105         dzoom = dzooma * pi/180.
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
112        xtild(i) = - pi + REAL(i) * depi /nmax2
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 ! ***  '
180        CALL ABORT_GCM("FXHYP", "", 1)
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
258      xlon2 = - pi + (REAL(i) + xuv - decalx) * depi / REAL(iim)
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
303       xxprim(i) = depi/ ( REAL(iim) * Xprimin )
304       xvrai(i)  =  xi + xzoom
305
3061500   CONTINUE
307
308
309
310       IF(ik.EQ.1.and.grossism.EQ.1.)  THEN
311         xvrai(1)    = xvrai(iip1)-depi
312         xxprim(1)   = xxprim(iip1)
313       ENDIF
314
315       DO i = 1 , iim
316        xlon(i)     = xvrai(i)
317        xprimm(i)   = xxprim(i)
318       ENDDO
319       DO i = 1, iim -1
320        IF( xvrai(i+1). LT. xvrai(i) )  THEN
321         WRITE(6,*) ' PBS. avec rlonu(',i+1,') plus petit que rlonu(',i,
322     ,  ')'
323        STOP 7
324        ENDIF
325       ENDDO
326c
327c   ... Reorganisation  des  longitudes  pour les avoir  entre - pi et pi ..
328c   ........................................................................
329
330       champmin =  1.e12
331       champmax = -1.e12
332       DO i = 1, iim
333        champmin = MIN( champmin,xvrai(i) )
334        champmax = MAX( champmax,xvrai(i) )
335       ENDDO
336
337      IF(champmin .GE.-pi-0.10.and.champmax.LE.pi+0.10 )  THEN
338                GO TO 1600
339      ELSE
340       WRITE(6,*) 'Reorganisation des longitudes pour avoir entre - pi',
341     ,  ' et pi '
342c
343        IF( xzoom.LE.0.)  THEN
344          IF( ik.EQ. 1 )  THEN
345          DO i = 1, iim
346           IF( xvrai(i).GE. - pi )  GO TO 80
347          ENDDO
348            WRITE(6,*)  ' PBS. 1 !  Xvrai plus petit que  - pi ! '
349            STOP 8
350 80       CONTINUE
351          is2 = i
352          ENDIF
353
354          IF( is2.NE. 1 )  THEN
355            DO ii = is2 , iim
356             xlon  (ii-is2+1) = xvrai(ii)
357             xprimm(ii-is2+1) = xxprim(ii)
358            ENDDO
359            DO ii = 1 , is2 -1
360             xlon  (ii+iim-is2+1) = xvrai(ii) + depi
361             xprimm(ii+iim-is2+1) = xxprim(ii)
362            ENDDO
363          ENDIF
364        ELSE
365          IF( ik.EQ.1 )  THEN
366           DO i = iim,1,-1
367             IF( xvrai(i).LE. pi ) GO TO 90
368           ENDDO
369             WRITE(6,*) ' PBS.  2 ! Xvrai plus grand  que   pi ! '
370              STOP 9
371 90        CONTINUE
372            is2 = i
373          ENDIF
374           idif = iim -is2
375           DO ii = 1, is2
376            xlon  (ii+idif) = xvrai(ii)
377            xprimm(ii+idif) = xxprim(ii)
378           ENDDO
379           DO ii = 1, idif
380            xlon (ii)  = xvrai (ii+is2) - depi
381            xprimm(ii) = xxprim(ii+is2)
382           ENDDO
383         ENDIF
384      ENDIF
385c
386c     .........   Fin  de la reorganisation   ............................
387
388 1600    CONTINUE
389
390
391         xlon  ( iip1)  = xlon(1) + depi
392         xprimm( iip1 ) = xprimm (1 )
393       
394         DO i = 1, iim+1
395         xvrai(i) = xlon(i)*180./pi
396         ENDDO
397
398         IF( ik.EQ.1 )  THEN
399c          WRITE(6,*)  ' XLON aux pts. V-0.25   apres ( en  deg. ) '
400c          WRITE(6,18)
401c          WRITE(6,68) xvrai
402c          WRITE(6,*) ' XPRIM k ',ik
403c          WRITE(6,566)  xprimm
404
405           DO i = 1,iim +1
406             rlonm025(i) = xlon( i )
407            xprimm025(i) = xprimm(i)
408           ENDDO
409         ELSE IF( ik.EQ.2 )  THEN
410c          WRITE(6,18)
411c          WRITE(6,*)  ' XLON aux pts. V   apres ( en  deg. ) '
412c          WRITE(6,68) xvrai
413c          WRITE(6,*) ' XPRIM k ',ik
414c          WRITE(6,566)  xprimm
415
416           DO i = 1,iim + 1
417             rlonv(i) = xlon( i )
418            xprimv(i) = xprimm(i)
419           ENDDO
420
421         ELSE IF( ik.EQ.3)   THEN
422c          WRITE(6,18)
423c          WRITE(6,*)  ' XLON aux pts. U   apres ( en  deg. ) '
424c          WRITE(6,68) xvrai
425c          WRITE(6,*) ' XPRIM ik ',ik
426c          WRITE(6,566)  xprimm
427
428           DO i = 1,iim + 1
429             rlonu(i) = xlon( i )
430            xprimu(i) = xprimm(i)
431           ENDDO
432
433         ELSE IF( ik.EQ.4 )  THEN
434c          WRITE(6,18)
435c          WRITE(6,*)  ' XLON aux pts. V+0.25   apres ( en  deg. ) '
436c          WRITE(6,68) xvrai
437c          WRITE(6,*) ' XPRIM ik ',ik
438c          WRITE(6,566)  xprimm
439
440           DO i = 1,iim + 1
441             rlonp025(i) = xlon( i )
442            xprimp025(i) = xprimm(i)
443           ENDDO
444
445         ENDIF
446
4475000    CONTINUE
448c
449       WRITE(6,18)
450c
451c    ...........  fin  de la boucle  do 5000      ............
452
453        DO i = 1, iim
454         xlon(i) = rlonv(i+1) - rlonv(i)
455        ENDDO
456        champmin =  1.e12
457        champmax = -1.e12
458        DO i = 1, iim
459         champmin = MIN( champmin, xlon(i) )
460         champmax = MAX( champmax, xlon(i) )
461        ENDDO
462         champmin = champmin * 180./pi
463         champmax = champmax * 180./pi
464
46518     FORMAT(/)
46624     FORMAT(2x,'Parametres xzoom,gross,tau ,dzoom pour fxhyp ',4f8.3)
46768     FORMAT(1x,7f9.2)
468566    FORMAT(1x,7f9.4)
469
470       RETURN
471       END
Note: See TracBrowser for help on using the repository browser.