source: LMDZ5/branches/LF-private/libf/dyn3dmem/fxhyp.F @ 5456

Last change on this file since 5456 was 1658, checked in by Laurent Fairhead, 12 years ago

Phasage de la dynamique parallele localisee (petite memoire) avec le tronc LMDZ4 (HEAD)
Validation effectuee par comparaison des fichiers de sorties debug (u, v, t, q, masse, etc ...) d'une simulation sans physique
faite avec la version du modele donnee par Y. Meurdesoif et la version phasee avec la r1428 (fin du tronc LMDZ4)


Phasing of the localised (low memory) parallel dynamics package with the LMDZ4 trunk version of LMDZ
Validation consisted in comparing output debug files (u, v, t, q, masse, etc... ) of a no physics simulation
run with the version of the code given by Y. Meurdesoif and this version phased with r1428 (HEAD of the LMDZ4 trunk)

File size: 11.7 KB
Line 
1!
2! $Id: fxhyp.F 1403 2010-07-01 09:02:53Z fairhead $
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*8 xlon(iip1),xprimm(iip1),xuv
51       REAL*8 xtild(0:nmax2)
52       REAL*8 fhyp(0:nmax2),ffdx,beta,Xprimt(0:nmax2)
53       REAL*8 Xf(0:nmax2),xxpr(0:nmax2)
54       REAL*8 xvrai(iip1),xxprim(iip1)
55       REAL*8 pi,depi,epsilon,xzoom,fa,fb
56       REAL*8 Xf1, Xfi , a0,a1,a2,a3,xi2
57       INTEGER i,it,ik,iter,ii,idif,ii1,ii2
58       REAL*8 xi,xo1,xmoy,xlon2,fxm,Xprimin
59       REAL*8 champmin,champmax,decalx
60       INTEGER is2
61       SAVE is2
62
63       REAL*8 heavyside
64
65       pi       = 2. * ASIN(1.)
66       depi     = 2. * pi
67       epsilon  = 1.e-3
68       xzoom    = xzoomdeg * pi/180.
69c
70           decalx   = .75
71       IF( grossism.EQ.1..AND.scal180 )  THEN
72           decalx   = 1.
73       ENDIF
74
75       WRITE(6,*) 'FXHYP scal180,decalx', scal180,decalx
76c
77       IF( dzooma.LT.1.)  THEN
78         dzoom = dzooma * depi
79       ELSEIF( dzooma.LT. 25. ) THEN
80         WRITE(6,*) ' Le param. dzoomx pour fxhyp est trop petit ! L aug
81     ,menter et relancer ! '
82         STOP 1
83       ELSE
84         dzoom = dzooma * pi/180.
85       ENDIF
86
87       WRITE(6,*) ' xzoom( rad.),grossism,tau,dzoom (radians)'
88       WRITE(6,24) xzoom,grossism,tau,dzoom
89
90       DO i = 0, nmax2
91        xtild(i) = - pi + REAL(i) * depi /nmax2
92       ENDDO
93
94       DO i = nmax, nmax2
95
96       fa  = tau*  ( dzoom/2.  - xtild(i) )
97       fb  = xtild(i) *  ( pi - xtild(i) )
98
99         IF( 200.* fb .LT. - fa )   THEN
100           fhyp ( i) = - 1.
101         ELSEIF( 200. * fb .LT. fa ) THEN
102           fhyp ( i) =   1.
103         ELSE
104            IF( ABS(fa).LT.1.e-13.AND.ABS(fb).LT.1.e-13)  THEN
105                IF(   200.*fb + fa.LT.1.e-10 )  THEN
106                    fhyp ( i ) = - 1.
107                ELSEIF( 200.*fb - fa.LT.1.e-10 )  THEN
108                    fhyp ( i )  =   1.
109                ENDIF
110            ELSE
111                    fhyp ( i )  =  TANH ( fa/fb )
112            ENDIF
113         ENDIF
114        IF ( xtild(i).EQ. 0. )  fhyp(i) =  1.
115        IF ( xtild(i).EQ. pi )  fhyp(i) = -1.
116
117       ENDDO
118
119cc  ....  Calcul  de  beta  ....
120
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
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
143       ENDIF
144
145       IF ( xmoy.EQ. 0. )  fxm  =  1.
146       IF ( xmoy.EQ. pi )  fxm  = -1.
147
148       ffdx = ffdx + fxm * ( xtild(i) - xtild(i-1) )
149
150       ENDDO
151
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
161c
162c   .....  calcul  de  Xprimt   .....
163c
164       
165       DO i = nmax, nmax2
166        Xprimt(i) = beta  + ( grossism - beta ) * fhyp(i)
167       ENDDO
168c   
169       DO i =  nmax+1, nmax2
170        Xprimt( nmax2 - i ) = Xprimt( i )
171       ENDDO
172c
173
174c   .....  Calcul  de  Xf     ........
175
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
196       ENDDO
197
198       DO i = nmax+1, nmax2
199        xxpr(nmax2-i+1) = xxpr(i)
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
213      WRITE(6,18)
214c
215      DO 5000  ik = 1, 4
216
217       IF( ik.EQ.1 )        THEN
218         xuv =  -0.25
219       ELSE IF ( ik.EQ.2 )  THEN
220         xuv =   0.
221       ELSE IF ( ik.EQ.3 )  THEN
222         xuv =   0.50
223       ELSE IF ( ik.EQ.4 )  THEN
224         xuv =   0.25
225       ENDIF
226
227      xo1   = 0.
228
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
235      DO 1500 i = ii1, ii2
236
237      xlon2 = - pi + (REAL(i) + xuv - decalx) * depi / REAL(iim)
238
239      Xfi    = xlon2
240c
241      DO 250 it =  nmax2,0,-1
242      IF( Xfi.GE.Xf(it))  GO TO 350
243250   CONTINUE
244
245      it = 0
246
247350   CONTINUE
248
249c    ......  Calcul de   Xf(xi)    ......
250c
251      xi  = xtild(it)
252
253      IF(it.EQ.nmax2)  THEN
254       it       = nmax2 -1
255       Xf(it+1) = pi
256      ENDIF
257c  .....................................................................
258c
259c   Appel de la routine qui calcule les coefficients a0,a1,a2,a3 d'un
260c   polynome de degre 3  qui passe  par les points (Xf(it),xtild(it) )
261c          et (Xf(it+1),xtild(it+1) )
262
263       CALL coefpoly ( Xf(it),Xf(it+1),Xprimt(it),Xprimt(it+1),
264     ,                xtild(it),xtild(it+1),  a0, a1, a2, a3  )
265
266       Xf1     = Xf(it)
267       Xprimin = a1 + 2.* a2 * xi + 3.*a3 * xi *xi
268
269       DO 500 iter = 1,300
270        xi = xi - ( Xf1 - Xfi )/ Xprimin
271
272        IF( ABS(xi-xo1).LE.epsilon)  GO TO 550
273         xo1      = xi
274         xi2      = xi * xi
275         Xf1      = a0 +  a1 * xi +     a2 * xi2  +     a3 * xi2 * xi
276         Xprimin  =       a1      + 2.* a2 *  xi  + 3.* a3 * xi2
277500   CONTINUE
278        WRITE(6,*) ' Pas de solution ***** ',i,xlon2,iter
279          STOP 6
280550   CONTINUE
281
282       xxprim(i) = depi/ ( REAL(iim) * Xprimin )
283       xvrai(i)  =  xi + xzoom
284
2851500   CONTINUE
286
287
288       IF(ik.EQ.1.and.grossism.EQ.1.)  THEN
289         xvrai(1)    = xvrai(iip1)-depi
290         xxprim(1)   = xxprim(iip1)
291       ENDIF
292       DO i = 1 , iim
293        xlon(i)     = xvrai(i)
294        xprimm(i)   = xxprim(i)
295       ENDDO
296       DO i = 1, iim -1
297        IF( xvrai(i+1). LT. xvrai(i) )  THEN
298         WRITE(6,*) ' PBS. avec rlonu(',i+1,') plus petit que rlonu(',i,
299     ,  ')'
300        STOP 7
301        ENDIF
302       ENDDO
303c
304c   ... Reorganisation  des  longitudes  pour les avoir  entre - pi et pi ..
305c   ........................................................................
306
307       champmin =  1.e12
308       champmax = -1.e12
309       DO i = 1, iim
310        champmin = MIN( champmin,xvrai(i) )
311        champmax = MAX( champmax,xvrai(i) )
312       ENDDO
313
314      IF(champmin .GE.-pi-0.10.and.champmax.LE.pi+0.10 )  THEN
315                GO TO 1600
316      ELSE
317       WRITE(6,*) 'Reorganisation des longitudes pour avoir entre - pi',
318     ,  ' et pi '
319c
320        IF( xzoom.LE.0.)  THEN
321          IF( ik.EQ. 1 )  THEN
322          DO i = 1, iim
323           IF( xvrai(i).GE. - pi )  GO TO 80
324          ENDDO
325            WRITE(6,*)  ' PBS. 1 !  Xvrai plus petit que  - pi ! '
326            STOP 8
327 80       CONTINUE
328          is2 = i
329          ENDIF
330
331          IF( is2.NE. 1 )  THEN
332            DO ii = is2 , iim
333             xlon  (ii-is2+1) = xvrai(ii)
334             xprimm(ii-is2+1) = xxprim(ii)
335            ENDDO
336            DO ii = 1 , is2 -1
337             xlon  (ii+iim-is2+1) = xvrai(ii) + depi
338             xprimm(ii+iim-is2+1) = xxprim(ii)
339            ENDDO
340          ENDIF
341        ELSE
342          IF( ik.EQ.1 )  THEN
343           DO i = iim,1,-1
344             IF( xvrai(i).LE. pi ) GO TO 90
345           ENDDO
346             WRITE(6,*) ' PBS.  2 ! Xvrai plus grand  que   pi ! '
347              STOP 9
348 90        CONTINUE
349            is2 = i
350          ENDIF
351           idif = iim -is2
352           DO ii = 1, is2
353            xlon  (ii+idif) = xvrai(ii)
354            xprimm(ii+idif) = xxprim(ii)
355           ENDDO
356           DO ii = 1, idif
357            xlon (ii)  = xvrai (ii+is2) - depi
358            xprimm(ii) = xxprim(ii+is2)
359           ENDDO
360         ENDIF
361      ENDIF
362c
363c     .........   Fin  de la reorganisation   ............................
364
365 1600    CONTINUE
366
367
368         xlon  ( iip1)  = xlon(1) + depi
369         xprimm( iip1 ) = xprimm (1 )
370       
371         DO i = 1, iim+1
372         xvrai(i) = xlon(i)*180./pi
373         ENDDO
374
375         IF( ik.EQ.1 )  THEN
376c          WRITE(6,*)  ' XLON aux pts. V-0.25   apres ( en  deg. ) '
377c          WRITE(6,18)
378c          WRITE(6,68) xvrai
379c          WRITE(6,*) ' XPRIM k ',ik
380c          WRITE(6,566)  xprimm
381
382           DO i = 1,iim +1
383             rlonm025(i) = xlon( i )
384            xprimm025(i) = xprimm(i)
385           ENDDO
386         ELSE IF( ik.EQ.2 )  THEN
387c          WRITE(6,18)
388c          WRITE(6,*)  ' XLON aux pts. V   apres ( en  deg. ) '
389c          WRITE(6,68) xvrai
390c          WRITE(6,*) ' XPRIM k ',ik
391c          WRITE(6,566)  xprimm
392
393           DO i = 1,iim + 1
394             rlonv(i) = xlon( i )
395            xprimv(i) = xprimm(i)
396           ENDDO
397
398         ELSE IF( ik.EQ.3)   THEN
399c          WRITE(6,18)
400c          WRITE(6,*)  ' XLON aux pts. U   apres ( en  deg. ) '
401c          WRITE(6,68) xvrai
402c          WRITE(6,*) ' XPRIM ik ',ik
403c          WRITE(6,566)  xprimm
404
405           DO i = 1,iim + 1
406             rlonu(i) = xlon( i )
407            xprimu(i) = xprimm(i)
408           ENDDO
409
410         ELSE IF( ik.EQ.4 )  THEN
411c          WRITE(6,18)
412c          WRITE(6,*)  ' XLON aux pts. V+0.25   apres ( en  deg. ) '
413c          WRITE(6,68) xvrai
414c          WRITE(6,*) ' XPRIM ik ',ik
415c          WRITE(6,566)  xprimm
416
417           DO i = 1,iim + 1
418             rlonp025(i) = xlon( i )
419            xprimp025(i) = xprimm(i)
420           ENDDO
421
422         ENDIF
423
4245000    CONTINUE
425c
426       WRITE(6,18)
427c
428c    ...........  fin  de la boucle  do 5000      ............
429
430        DO i = 1, iim
431         xlon(i) = rlonv(i+1) - rlonv(i)
432        ENDDO
433        champmin =  1.e12
434        champmax = -1.e12
435        DO i = 1, iim
436         champmin = MIN( champmin, xlon(i) )
437         champmax = MAX( champmax, xlon(i) )
438        ENDDO
439         champmin = champmin * 180./pi
440         champmax = champmax * 180./pi
441
44218     FORMAT(/)
44324     FORMAT(2x,'Parametres xzoom,gross,tau ,dzoom pour fxhyp ',4f8.3)
44468     FORMAT(1x,7f9.2)
445566    FORMAT(1x,7f9.4)
446
447       RETURN
448       END
Note: See TracBrowser for help on using the repository browser.