source: LMDZ4/trunk/libf/dyn3dpar/fxhyp.F

Last change on this file was 1403, checked in by Laurent Fairhead, 14 years ago

Merged LMDZ4V5.0-dev branch changes r1292:r1399 to trunk.

Validation:
Validation consisted in compiling the HEAD revision of the trunk,
LMDZ4V5.0-dev branch and the merged sources and running different
configurations on local and SX8 machines comparing results.

Local machine: bench configuration, 32x24x11, gfortran

  • IPSLCM5A configuration (comparison between trunk and merged sources):
    • numerical convergence on dynamical fields over 3 days
    • start files are equivalent (except for RN and PB fields)
    • daily history files equivalent
  • MH07 configuration, new physics package (comparison between LMDZ4V5.0-dev branch and merged sources):
    • numerical convergence on dynamical fields over 3 days
    • start files are equivalent (except for RN and PB fields)
    • daily history files equivalent

SX8 machine (brodie), 96x95x39 on 4 processors:

  • IPSLCM5A configuration:
    • start files are equivalent (except for RN and PB fields)
    • monthly history files equivalent
  • MH07 configuration:
    • start files are equivalent (except for RN and PB fields)
    • monthly history files equivalent

Changes to the makegcm and create_make_gcm scripts to take into account
main programs in F90 files

Fusion de la branche LMDZ4V5.0-dev (r1292:r1399) au tronc principal

Validation:
La validation a consisté à compiler la HEAD de le trunk et de la banche
LMDZ4V5.0-dev et les sources fusionnées et de faire tourner le modéle selon
différentes configurations en local et sur SX8 et de comparer les résultats

En local: 32x24x11, config bench/gfortran

  • pour une config IPSLCM5A (comparaison tronc/fusion):
    • convergence numérique sur les champs dynamiques après 3 jours
    • restart et restartphy égaux (à part sur RN et Pb)
    • fichiers histoire égaux
  • pour une config nlle physique (MH07) (comparaison LMDZ4v5.0-dev/fusion):
    • convergence numérique sur les champs dynamiques après 3 jours
    • restart et restartphy égaux
    • fichiers histoire équivalents

Sur brodie, 96x95x39 sur 4 proc:

  • pour une config IPSLCM5A:
    • restart et restartphy égaux (à part sur RN et PB)
    • pas de différence dans les fichiers histmth.nc
  • pour une config MH07
    • restart et restartphy égaux (à part sur RN et PB)
    • pas de différence dans les fichiers histmth.nc

Changement sur makegcm et create_make-gcm pour pouvoir prendre en compte des
programmes principaux en *F90
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.7 KB
Line 
1!
2! $Id: fxhyp.F 1403 2010-07-01 09:02:53Z aclsce $
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.