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

Last change on this file since 204 was 203, checked in by lmdz, 24 years ago

Mise en oeuvre du nouveau zoom. Modifs P. 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: 10.8 KB
Line 
1c
2c $Header
3c
4       SUBROUTINE fxhyp ( xzoomdeg,grossism,dzoom,tau ,
5     , rlonm025,xprimm025,rlonv,xprimv,rlonu,xprimu,rlonp025,xprimp025)
6
7c      Auteur :  P. Le Van
8
9       IMPLICIT NONE
10
11c    Calcule les longitudes et derivees dans la grille du GCM pour une
12c     fonction f(x) a tangente  hyperbolique  .
13c
14c     grossism etant le grossissement ( = 2 si 2 fois, = 3 si 3 fois,etc.)
15c     dzoom  etant  la distance totale de la zone du zoom
16c     tau  la raideur de la transition de l'interieur a l'exterieur du zoom
17c
18c    On doit avoir grossism x dzoom <  pi ( radians )   , en longitude.
19c   ********************************************************************
20
21
22       INTEGER nmax, nmax2
23       PARAMETER (  nmax = 30000, nmax2 = 2*nmax )
24
25#include "dimensions.h"
26#include "paramet.h"
27       
28c     ......  arguments  d'entree   .......
29c
30       REAL xzoomdeg,dzoom,tau,grossism
31
32c    ......   arguments  de  sortie  ......
33
34       REAL rlonm025(iip1),xprimm025(iip1),rlonv(iip1),xprimv(iip1),
35     ,  rlonu(iip1),xprimu(iip1),rlonp025(iip1),xprimp025(iip1)
36
37c     .... variables locales  ....
38c
39       REAL*8 xlon(iip1),xprimm(iip1),xuv
40       REAL*8 xtild(0:nmax2)
41       REAL*8 fhyp(0:nmax),ffdx,beta,Xprimt(0:nmax2)
42       REAL*8 Xf(0:nmax2),xxpr(0:nmax2)
43       REAL*8 xvrai(iip1),xxprim(iip1)
44       REAL*8 pi,depi,epsilon,xzoom,fa,fb
45       REAL*8 Xf1, Xfi , a0,a1,a2,a3,xi2
46       INTEGER i,it,ik,iter,ii,idif
47       REAL*8 xi,xo1,xmoy,xlon2,fxm,Xprimin
48       REAL*8 champmin,champmax
49       INTEGER is2
50       SAVE is2
51       REAL*8 heavyside
52       EXTERNAL coefpoly,heavyside
53
54       pi       = 2. * ASIN(1.)
55       depi     = 2. * pi
56       epsilon  = 1.e-3
57       xzoom    = xzoomdeg * pi/180.
58
59       WRITE(6,24) xzoomdeg,grossism,tau,dzoom
60       IF( dzoom.LT.1.)  THEN
61         dzoom = dzoom * depi
62       ELSEIF( dzoom.LT. 25. ) THEN
63         WRITE(6,*) ' Le param. dzoomy pour fxhyp est trop petit ! L aug
64     ,menter et relancer ! '
65         STOP 1
66       ELSE
67         dzoom = dzoom * pi/180.
68       ENDIF
69
70       DO i = 0, nmax2
71        xtild(i) = - pi + FLOAT(i) * depi /nmax2
72       ENDDO
73
74       DO i = 0, nmax2
75
76       fa  = tau*  ( dzoom/2.  - xtild(i) )
77       fb  = xtild(i) *  ( pi - xtild(i) )
78
79
80       IF( 200.* fb .LT. - fa )   THEN
81         fhyp ( i) = - 1.
82       ELSEIF( 200. * fb .LT. fa ) THEN
83         fhyp ( i) =   1.
84       ELSE
85         fhyp(i) =  TANH ( fa/fb )
86       ENDIF
87
88       IF ( xtild(i).EQ. 0. )  fhyp(i) =  1.
89       IF ( xtild(i).EQ. pi )  fhyp(i) = -1.
90
91       ENDDO
92
93cc  ....  Calcul  de  beta  ....
94c   ............................
95
96       ffdx = 0.
97
98       DO i = nmax +1,nmax2
99
100       xmoy    = 0.5 * ( xtild(i-1) + xtild( i ) )
101       fa  = tau*  ( dzoom/2.  - xmoy )
102       fb  = xmoy *  ( pi - xmoy )
103
104       IF( 200.* fb .LT. - fa )   THEN
105         fxm = - 1.
106       ELSEIF( 200. * fb .LT. fa ) THEN
107         fxm =   1.
108       ELSE
109         fxm =  TANH ( fa/fb )
110       ENDIF
111
112       IF ( xmoy.EQ. 0. )  fhyp(i) =  1.
113       IF ( xmoy.EQ. pi )  fhyp(i) = -1.
114       ffdx = ffdx + fxm * ( xtild(i) - xtild(i-1) )
115
116       ENDDO
117
118        beta  = ( grossism * ffdx - pi ) / ( ffdx - pi )
119ccc        WRITE(6,*) ' ** X  beta **',beta,grossism
120
121       IF( 2.*beta - grossism.LE. 0.)  THEN
122        WRITE(6,*) ' **  Attention ! La valeur beta calculee dans la rou
123     ,tine fxhyp est mauvaise ! '
124        WRITE(6,*)'Modifier les valeurs de  grossismx ,tau ou dzoomx ',
125     , ' et relancer ! ***  '
126        CALL ABORT
127       ENDIF
128c
129c   .....  calcul  de  Xprimt   .....
130c
131       
132       DO i = nmax, nmax2
133        Xprimt(i) = beta  + ( grossism - beta ) * fhyp(i)
134       ENDDO
135c   
136       DO i =  nmax+1, nmax2
137        Xprimt( nmax2 - i ) = Xprimt( i )
138       ENDDO
139c
140
141c   .....  Calcul  de  Xf     ........
142
143       Xf(0) = - pi
144
145       DO i =  nmax +1, nmax2
146
147       xmoy    = 0.5 * ( xtild(i-1) + xtild( i ) )
148       fa  = tau*  ( dzoom/2.  - xmoy )
149       fb  = xmoy *  ( pi - xmoy )
150
151       IF( 200.* fb .LT. - fa )   THEN
152         fxm = - 1.
153       ELSEIF( 200. * fb .LT. fa ) THEN
154         fxm =   1.
155       ELSE
156         fxm =  TANH ( fa/fb )
157       ENDIF
158
159       IF ( xmoy.EQ. 0. )  fxm =  1.
160       IF ( xmoy.EQ. pi )  fxm = -1.
161       xxpr(i)    = beta + ( grossism - beta ) * fxm
162
163       ENDDO
164
165       DO i = nmax+1, nmax2
166        xxpr(nmax2-i+1) = xxpr(i)
167       ENDDO
168
169        DO i=1,nmax2
170         Xf(i)   = Xf(i-1) + xxpr(i) * ( xtild(i) - xtild(i-1) )
171        ENDDO
172
173
174c    *****************************************************************
175c
176
177c     .....  xuv = 0.   si  calcul  aux pts   scalaires   ........
178c     .....  xuv = 0.5  si  calcul  aux pts      U        ........
179c
180c
181      DO 5000  ik = 1, 4
182
183       IF( ik.EQ.1 )        THEN
184         xuv = - 0.25
185       ELSE IF ( ik.EQ.2 )  THEN
186         xuv =   0.
187       ELSE IF ( ik.EQ.3 )  THEN
188         xuv =   0.5
189       ELSE IF ( ik.EQ.4 )  THEN
190         xuv =   0.25
191       ENDIF
192
193      xo1   = 0.
194
195      DO 1500 i = 1, iim
196
197      xlon2 = - pi + (FLOAT(i) + xuv - 0.75) * depi / FLOAT(iim)
198
199      Xfi    = xlon2
200c
201      DO 250 it =  nmax2,0,-1
202      IF( Xfi.GE.Xf(it))  GO TO 350
203250   CONTINUE
204
205      it = 0
206
207350   CONTINUE
208
209c    ......  Calcul de   Xf(xi)    ......
210c
211      xi  = xtild(it)
212
213      IF(it.EQ.nmax2)  THEN
214       it       = nmax2 -1
215       Xf(it+1) = pi
216      ENDIF
217c  .....................................................................
218c
219c   Appel de la routine qui calcule les coefficients a0,a1,a2,a3 d'un
220c   polynome de degre 3  qui passe  par les points (Xf(it),xtild(it) )
221c          et (Xf(it+1),xtild(it+1) )
222
223       CALL coefpoly ( Xf(it),Xf(it+1),Xprimt(it),Xprimt(it+1),
224     ,                xtild(it),xtild(it+1),  a0, a1, a2, a3  )
225
226       Xf1     = Xf(it)
227       Xprimin = a1 + 2.* a2 * xi + 3.*a3 * xi *xi
228
229       DO 500 iter = 1,300
230        xi = xi - ( Xf1 - Xfi )/ Xprimin
231
232        IF( ABS(xi-xo1).LE.epsilon)  GO TO 550
233         xo1      = xi
234         xi2      = xi * xi
235         Xf1      = a0 +  a1 * xi +     a2 * xi2  +     a3 * xi2 * xi
236         Xprimin  =       a1      + 2.* a2 *  xi  + 3.* a3 * xi2
237500   CONTINUE
238        WRITE(6,*) ' Pas de solution ***** ',i,xlon2,iter
239          STOP 6
240550   CONTINUE
241
242       xxprim(i) = depi/ ( FLOAT(iim) * Xprimin )
243       xvrai(i)  =  xi + xzoom
244
2451500   CONTINUE
246
247       DO i = 1 , iim
248        xlon(i)     = xvrai(i)
249        xprimm(i)   = xxprim(i)
250       ENDDO
251       
252       DO i = 1, iim -1
253        IF( xvrai(i+1). LT. xvrai(i) )  THEN
254         WRITE(6,*) ' PBS. avec rlonu(',i+1,') plus petit que rlonu(',i,
255     ,  ')'
256        STOP 7
257        ENDIF
258       ENDDO
259c
260c   ... Reorganisation  des  longitudes  pour les avoir  entre - pi et pi ..
261c   ........................................................................
262
263       champmin =  1.e12
264       champmax = -1.e12
265       DO i = 1, iim
266        champmin = MIN( champmin,xvrai(i) )
267        champmax = MAX( champmax,xvrai(i) )
268       ENDDO
269
270      IF(champmin .GE. - pi .AND. champmax.LE. pi )  THEN
271                GO TO 1600
272      ELSE
273       WRITE(6,18)
274       WRITE(6,*) 'Reorganisation des longitudes pour avoir entre - pi',
275     ,  ' et pi '
276c
277        IF( xzoom.LE.0.)  THEN
278          IF( ik.EQ. 1 )  THEN
279          DO i = 1, iim
280           IF( xvrai(i).GE. - pi )  GO TO 80
281          ENDDO
282            WRITE(6,*)  ' PBS. 1 !  Xvrai plus petit que  - pi ! '
283            STOP 8
284 80       CONTINUE
285          is2 = i
286          ENDIF
287
288          IF( is2.NE. 1 )  THEN
289            DO ii = is2 , iim
290             xlon  (ii-is2+1) = xvrai(ii)
291             xprimm(ii-is2+1) = xxprim(ii)
292            ENDDO
293            DO ii = 1 , is2 -1
294             xlon  (ii+iim-is2+1) = xvrai(ii) + depi
295             xprimm(ii+iim-is2+1) = xxprim(ii)
296            ENDDO
297          ENDIF
298        ELSE
299          IF( ik.EQ.1 )  THEN
300           DO i = iim,1,-1
301            IF( xvrai(i).LE. pi )  GO TO 90
302           ENDDO
303             WRITE(6,*) ' PBS.  2 ! Xvrai plus grand  que   pi ! '
304              STOP 9
305 90        CONTINUE
306            is2 = i
307          ENDIF
308           idif = iim -is2
309           DO ii = 1, is2
310            xlon  (ii+idif) = xvrai(ii)
311            xprimm(ii+idif) = xxprim(ii)
312           ENDDO
313           DO ii = 1, idif
314            xlon (ii)  = xvrai (ii+is2) - depi
315            xprimm(ii) = xxprim(ii+is2)
316           ENDDO
317         ENDIF
318      ENDIF
319c
320c     .........   Fin  de la reorganisation   ............................
321
322 1600    CONTINUE
323
324
325         xlon  ( iip1)  = xlon(1) + depi
326         xprimm( iip1 ) = xprimm (1 )
327       
328         DO i = 1, iim+1
329         xvrai(i) = xlon(i)*180./pi
330         ENDDO
331
332         IF( ik.EQ.1 )  THEN
333          WRITE(6,*)  ' XLON aux pts. V-0.25   apres ( en  deg. ) '
334          WRITE(6,18)
335          WRITE(6,68) xvrai
336ccc          WRITE(6,*) ' XPRIM k ',ik
337ccc          WRITE(6,566)  xprimm
338
339           DO i = 1,iim + 1
340             rlonm025(i) = xlon( i )
341            xprimm025(i) = xprimm(i)
342           ENDDO
343
344         ELSE IF( ik.EQ.2 )  THEN
345          WRITE(6,18)
346          WRITE(6,*)  ' XLON aux pts. V   apres ( en  deg. ) '
347          WRITE(6,68) xvrai
348ccc          WRITE(6,*) ' XPRIM k ',ik
349ccc          WRITE(6,566)  xprimm
350
351           DO i = 1,iim + 1
352             rlonv(i) = xlon( i )
353            xprimv(i) = xprimm(i)
354           ENDDO
355
356         ELSE IF( ik.EQ.3)   THEN
357          WRITE(6,18)
358          WRITE(6,*)  ' XLON aux pts. U   apres ( en  deg. ) '
359          WRITE(6,68) xvrai
360ccc          WRITE(6,*) ' XPRIM ik ',ik
361ccc          WRITE(6,566)  xprimm
362
363           DO i = 1,iim + 1
364             rlonu(i) = xlon( i )
365            xprimu(i) = xprimm(i)
366           ENDDO
367
368         ELSE IF( ik.EQ.4 )  THEN
369          WRITE(6,18)
370          WRITE(6,*)  ' XLON aux pts. V+0.25   apres ( en  deg. ) '
371          WRITE(6,68) xvrai
372ccc          WRITE(6,*) ' XPRIM ik ',ik
373ccc          WRITE(6,566)  xprimm
374
375           DO i = 1,iim + 1
376             rlonp025(i) = xlon( i )
377            xprimp025(i) = xprimm(i)
378           ENDDO
379
380         ENDIF
381
3825000    CONTINUE
383c
384c    ...........  fin  de la boucle  do 5000      ............
385
386        DO i = 1, iim
387         xlon(i) = rlonv(i+1) - rlonv(i)
388        ENDDO
389        champmin =  1.e12
390        champmax = -1.e12
391        DO i = 1, iim
392         champmin = MIN( champmin, xlon(i) )
393         champmax = MAX( champmax, xlon(i) )
394        ENDDO
395         champmin = champmin * 180./pi
396         champmax = champmax * 180./pi
397
398        WRITE(6,18)
399        WRITE(6,*) '  Longitudes '
400        WRITE(6,18)
401        WRITE(6,3)  champmin, champmax
402        WRITE(6,*) ' Si cette derniere est trop lache , modifiez les par
403     ,ametres  grossism , tau , dzoom pour X et repasser ! '
404        WRITE(6,18)
405
4063      Format(1x, ' Au centre du zoom , la longueur de la maille est',
407     ,  ' d environ ',f8.2 ,' degres  ',
408     , ' alors que la maille en dehors de la zone du zoom est d environ
409     , ', f8.2,' degres ' )
41018     FORMAT(/)
41124     FORMAT(2x,' Parametres yzoom,gross,tau ,dzoom pour fyhyp ',4f8.3)
41268     FORMAT(1x,7f9.2)
413
414       RETURN
415       END
Note: See TracBrowser for help on using the repository browser.